Arctic sea ice, Eurasia snow, and extreme winter haze in China

See allHide authors and affiliations

Science Advances  15 Mar 2017:
Vol. 3, no. 3, e1602751
DOI: 10.1126/sciadv.1602751


The East China Plains (ECP) region experienced the worst haze pollution on record for January in 2013. We show that the unprecedented haze event is due to the extremely poor ventilation conditions, which had not been seen in the preceding three decades. Statistical analysis suggests that the extremely poor ventilation conditions are linked to Arctic sea ice loss in the preceding autumn and extensive boreal snowfall in the earlier winter. We identify the regional circulation mode that leads to extremely poor ventilation over the ECP region. Climate model simulations indicate that boreal cryospheric forcing enhances the regional circulation mode of poor ventilation in the ECP region and provides conducive conditions for extreme haze such as that of 2013. Consequently, extreme haze events in winter will likely occur at a higher frequency in China as a result of the changing boreal cryosphere, posing difficult challenges for winter haze mitigation but providing a strong incentive for greenhouse gas emission reduction.

  • extreme haze
  • pollution in China
  • boreal cryosphere changes
  • climate change and impacts


A consequence of rapid economic growth in China is a deterioration of air quality. An air pollution record was set in January 2013, with unprecedented large-scale haze lasting almost an entire month. During this so-called “airpocalypse” period, ~70% of the 74 major cities exceeded the daily PM2.5 (particulate matter with a diameter of 2.5 μm or less; see table S1 for a list of all acronyms) ambient air quality standard of China (75 μg m−3), with the maximum daily PM2.5 reaching 766 μg m−3 and the monthly mean concentration reaching as high as 130 μg m−3 (1). Exposure to such high PM concentrations endangers public health with increasing risks of cardiovascular and respiratory morbidity (2, 3).

Because no sudden rapid emission surge of natural or anthropogenic emissions in that month over eastern China was reported, stagnant meteorological conditions favoring the high aerosol formation (4) and accumulation (5, 6) appeared to be a major factor contributing to the extreme haze. Although recent studies indicated a decadal weakening trend of the East Asian winter monsoon (EAWM) (7) and consequently decreasing wind speed (8) and increasing aerosol concentrations (9), the underlying climate factors modulating short-term weather changes relevant to regional air quality are still not understood.

This study aims to place the occurrence of the January 2013 extreme haze in the context of historical ventilation conditions in the last 35 years. The region of East China Plains (ECPs; 112°E to 122°E, 30°N to 41°N; Fig. 1A) is the focus of this study. The region hosts a large portion of the Chinese population and suffers from severe air pollution problems. It resembles a horseshoe-shaped basin, where the ventilation of air pollutants relies on large-scale weather systems (10). Pollutant ventilation can be either horizontal or vertical. Using the NCEP (National Centers for Environmental Prediction)/NCAR (National Center for Atmospheric Research) reanalysis data (11), we first compute a normalized near-surface wind speed index (WSI) for horizontal ventilation and a potential air temperature gradient index (ATGI) for vertical ventilation and then construct a synthetic meteorological index—pollution potential index (PPI)—to better quantify the synergistic effect of ventilation on regional air pollution in January. We then explore the relationship between PPI and regional climate variability through comprehensive statistical approaches including maximum covariance analysis (MCA) (12) and principal components analysis (PCA) (12). Using these analyses as a guide, we investigate how anomalous synoptic patterns under the background of changing climate contribute to poor ventilation and extreme haze in China. We also conduct sensitivity simulations using the state-of-the-science Community Earth System Model (CESM) and examine the modeling results from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) (13) to validate the statistical findings.

Fig. 1 PM pollution and ventilation conditions over East Asia.

(A) 2013 monthly Moderate Resolution Imaging Spectroradiometer (MODIS) AOD (unitless) at 550 nm onboard Aqua satellite; (B) time series of aerosol observations, PPI, and CFI; their correlation coefficients (r values) with PPI are shown in parentheses; (C) 2013 distributions of normalized surface WSI (unitless); (D) 2013 distributions of normalized potential ATGI (unitless). In (C) and (D), black dots (crosses) denote the 99% (95%) significance level based on the bootstrapping method. The red rectangular box in (A) and the black box in (C) and (D) show the ECP region. All results are for January.


Poor ventilation and extreme haze in China

In January 2013, the ECP region was characterized by large negative surface wind speed anomalies (weakened horizontal dispersion; Fig. 1C) and positive potential air temperature gradient anomalies (more stable atmosphere with weakened vertical convection; Fig. 1D). Although precipitation is also an important factor to remove PM through wet scavenging, previous studies indicated less significant dependence of PM concentrations on precipitation in winter (14), which was corroborated in our analysis. As expected, the regionally averaged ventilation indices (WSI and ATGI) are correlated with historical haze observations including PM10 over the ECP region, PM2.5 in Beijing, visibility inverse (ViI) (15), and satellite column aerosol optical depth (AOD) data (table S2) (16). To simplify the multivariate statistical analysis, we construct a synthetic meteorological index, PPI, as a correlation weighted average of WSI and ATGI (hereafter, WSI, ATGI, and PPI are regionally averaged over the ECP region if not noted otherwise; see Materials and Methods for details). The PM10 data suggest roughly equal weighting of horizontal ventilation (WSI) and vertical ventilation (ATGI) indices in PPI. The new proxy (PPI) correlates better with ground observations such as PM10 (r = 0.92), PM2.5 (r = 0.79), and ViI (r = 0.62) than satellite AOD data (r = 0.43 to 0.50; Fig. 1B and table S2), indicating that PPI is more representative for near-surface air quality than column aerosol loading. In a historical perspective, both calculated PPI and observed ViI show that the January 2013 extreme haze is unprecedented in the last 30 years (Fig. 1B).

The role of climate variability

To investigate the association of climate variability to the PPI extreme, we apply the PCA (12) to PPI, multiple atmospheric variables including the Arctic Oscillation index (17) and EAWM indices (1820), and three forcing factors, Arctic sea ice concentration (SIC; Fig. 1A) (21) in the preceding autumn, boreal Eurasia snow cover extent (SCE; Fig. 1A) (22) in early winter, and El Niño/Southern Oscillation (ENSO) (23), for the past three decades (1980 to 2015; see Materials and Methods and table S3 for details). The long-term records of PPI and climate variables are necessary for climate-related analysis. We note that the construction of PPI is not directly related to any atmospheric variable or forcing factor used here. Unlike the other years, all related PCs contribute positively to the 2013 extreme, and the large contributors PC5 and PC6 are correlated with SCE and SIC, respectively (fig. S1 and table S4). As in the case of PPI, we construct a normalized cryospheric forcing index (CFI), which is an average of correlation weighted SIC and SCE (see Materials and Methods for details). As such, we have one representative variable for cryospheric forcing discussion. CFI can explain 42% of the total PPI variance (table S2). In 2013, the large CFI extreme corresponds well to the PPI extreme (Fig. 1B).

To examine the atmospheric processes linking cryospheric changes to PPI, we apply the MCA (12), which is totally independent from the PCA, to analyze the association of mid-latitude to high-latitude synoptic weather patterns [the 850-hPa geopotential height (Z850) field of the reanalysis data] with PPI (see Materials and Methods for details). We note that the formulation of PPI is not directly related to the Z850 field, that the region of interest for the gridded PPI field (Fig. 2C) is of a smaller domain than that for the Z850 field (Fig. 2B), and that the MCA modes remain the same if the 2013 data are excluded from the analysis. The most dominant coupling modes of the Z850 and PPI fields account for >30% of the total covariance. The highest Z850 and PPI mode intensity also corresponds to the extreme PPI value in 2013 (Fig. 2A). The spatial pattern of the first MCA Z850 mode (Fig. 2B) resembles the circulation anomaly in 2013 (Fig. 2D). Similarly, the spatial pattern of the first MCA PPI mode resembles the PPI distribution in 2013 (Fig. 2C). Furthermore, the time series of the first MCA PPI mode intensity is highly correlated with the average PPI over the ECP region (r = 0.93). Therefore, the extreme condition of 2013 can be understood using the more general MCA modes; that is, the poor ventilation condition over the ECP region represented by the first MCA PPI mode is driven by the regional circulation pattern represented by the first MCA Z850 mode. In contrast to the climatology characterized by large pressure gradients between the continent and the oceans (Fig. 2B), the first MCA Z850 mode shows a reversed northeast-southwest pressure gradient with anticyclonic anomalies in the Arctic and northeast Asia and a cyclonic anomaly over central Siberia, leading to weakened monsoon wind and enhanced PPI over the ECP region.

Fig. 2 Influence of the regional circulation on PPI.

(A) The first MCA mode intensity (circles; unitless) of January PPI and Z850; color shading (unitless) denotes PPI values from 1981 to 2015; (B) the spatial pattern of the first Z850 MCA mode (color shading; unitless) and Z850 climatology (contour lines; in meters); (C) the spatial pattern of the first PPI MCA mode (color shading; unitless) and PPI fields (contour lines; unitless) in 2013; (D) the 2013 Z850 anomalies (color shading; in meters) and Z850 climatology (contour lines; in meters). In (D), black dots denote the 95% significance level based on the bootstrapping method; green rectangles denote the regions of subplots (B) and (C); H and L indicate the location of the Siberian High and the Aleutian Low, respectively. All results are for January.

The 2013 event is therefore a manifestation of the first MCA mode characterized by poor ventilation over the ECP region. It is plausible that cryospheric forcing due to Arctic sea ice and Eurasian snowfall identified in the PCA enhances the first MCA Z850 mode, leading to high PPI and hence heavy haze over the ECP region. To test this hypothesis, we designed sensitivity simulations using the state-of-the-science CESM (version 1.2.1). We conducted a 30-year control (CTRL) run with prescribed climatological Arctic sea ice and sea surface temperature (SST) (21) and three sensitivity experiments with 30 ensemble members each: the first sensitivity simulation set (SENS1) with climatological SIC and SST data replaced by Arctic SIC and SST data in August to November as observed in 2012; the second sensitivity simulation set (SENS2) with prescribed Eurasia snow water equivalent (SWE) data (24) in October and November as observed in 2012; and the third sensitivity simulation set (SENS3) with the observed SIC, SST, and SWE data, in the same manner as in SENS1 and SENS2 (see Materials and Methods and fig. S2).

The response of the PPI cumulative distribution function (CDF) to cryospheric forcing in the model is compared to the reanalysis data (Fig. 3). Relative to the CTRL simulations and reanalysis data, Arctic sea ice and Eurasian snow forcing clearly shift PPI toward positive values (Fig. 3A). The fraction of positive PPI increases from 49% in reanalysis data and 52% in CTRL simulations to 58% in SENS1 with sea ice forcing, 62% in SENS2 with snow forcing, and 66% in SENS3 with both sea ice and snow forcing. Consequently, the ensemble averaged value of PPI in winter increases from 0 in CTRL to 0.06 (P = 0.34) in SENS1 and 0.20 in both SENS2 (P = 0.005) and SENS3 (P = 0.003). To investigate the extreme cases more relevant to the 2013 event, we examined the distribution of sensitivity simulation extreme ensemble members, defined as the PPI values of which greater than the 95th percentile value of the CTRL ensemble (Embedded Image). The number of SENS1 extreme members is similar to that of CTRL, whereas the number of SENS2 and SENS3 extreme members increases by 80 and 100%, respectively. Eurasian snow forcing (in SENS2 and SENS3) appears to have a larger effect on the extreme PPI values, particularly over the ECP region, than Arctic sea ice forcing (SENS1) (Fig. 3, B to D), in agreement with the observations (table S2).

Fig. 3 PPI responses to cryosphere forcing in the CESM sensitivity simulations.

(A) The CDF of PPI over the ECP region in the reanalysis data and CESM sensitivity simulations. The red solid vertical line is at the PPI value of 0; (B) the spatial distribution of ensemble averaged PPI fields of the extreme members (Embedded Image) in SENS1; (C) same as (B) but in SENS2; (D) same as (B) but in SENS3. In (B) to (D), black dots denote the 99% significance level based on the bootstrapping method.

The intensity distributions of the first MCA modes of Z850 and PPI (Fig. 2) in the ensemble simulations are symmetrical in the CTRL run because of internal model variability (fig. S3). These modes clearly shift to higher intensities, toward the 2013 anomaly when Arctic sea ice and Eurasian snow forcing are added, confirming that the cryospheric forcing provides conducive conditions for extreme haze over the ECP region. These modes respond more to Arctic sea ice forcing than to Eurasian snow forcing (fig. S3), likely reflecting in part the difference between climate model simulations and reanalysis data.

Future projections from CMIP5

Observation-based PCA and MCA and CESM simulations in our study indicate the contributions of decreasing Arctic sea ice (25, 26) and increasing Eurasian winter snowfall (27) to poor ventilation over the ECP region. As such, the capability of climate models is challenged in reproducing historical observations and projecting future changes of winter ventilation and air pollution over eastern China. We analyzed the modeling results from the CMIP5 project (13) using CESM1-CAM5 (the Community Atmosphere Model version 5) and CCSM4 (the Community Climate System Model version 4) models. Although the SIC hindcast agrees well with the observations (Fig. 4A), the increase of cryospheric forcing is underestimated in the climate models (Fig. 4C), which is largely attributed to the model’s failure to reproduce the observed increase of early winter boreal snow (Fig. 4B). The maximum simulated correlation coefficient between simulated PPI and CFI barely reaches that between observed PPI and SIC (|r| = 0.43; Fig. 4E), which is much lower than the observed value of 0.65 in the last 35 years. The reason is that both the observations and CESM sensitivity simulations suggest that increasing snowfalls tend to increase PPI (Fig. 3 and table S2). However, the CMIP5-simulated early winter SCE over Eurasia shows a consistently decreasing trend (Fig. 4B), in contrast to the observed increasing trend (27), leading to a failure of the CMIP5 models to project high PPI values in response to the rapid increases of SCE and CFI (Fig. 4D). Consequently, these climate projections underestimate the climate penalty and the benefit of climate change mitigation on wintertime air quality over eastern China.

Fig. 4 Time series of cryospheric forcing factors, PPI, and the correlations of PPI with cryospheric forcing based on the observations, reanalysis data, and CMIP5 projections.

(A) Comparisons of SIC observations and ensemble averaged CESM1(CAM5) (three ensemble members) and CCSM4 (six ensemble members) CMIP5 simulations (unitless; color shading denotes ±1 SD); (B) same as (A) but for SCE; (C) same as (A) but for CFI; (D) same as (A) but for PPI; (E) correlation coefficients of PPI with SCE, SIC, and CFI based on the 35-year observations and reanalysis data [symbols; unitless; the absolute value of the PPI-SIC correlation coefficient (magenta square) is shown for comparison purpose] and 35-year moving correlation coefficients between simulated PPI and CFI of CMIP5 ensembles (lines; unitless; the data are plotted at the last year of each 35-year period).


The winter cold extremes and heavy snowfalls in northern high latitudes are partially attributable to Arctic sea ice loss (2830). In 2012, most boreal regions experienced a chilly early winter with anomalously heavy snowfalls after a record-breaking decline of Arctic sea ice in September 2012. The surface temperature anomalies are particularly apparent in northern Eurasia, with strong warming over the Arctic and cooling over boreal Eurasia in the subsequent months, and model sensitivity results suggest that its formation is attributable to both Arctic sea ice and Eurasia snow forcing (fig. S4). Previous studies suggest that stratosphere-troposphere coupling is critical in climate responses over mid-latitude regions to polar forcing (31, 32). Arctic sea ice and Eurasia snow forcing strengthen the vertical propagation of Rossby waves from the troposphere to the stratosphere and weaken the stratospheric polar vortex, leading to a lower tropospheric response resembling our dominant MCA mode of the geopotential height field (Fig. 2B) (32, 33). The CESM simulations (figs. S4 and S5) are also consistent with this mechanism.

If the decrease of Arctic sea ice and the increase of Eurasian snowfall continue, we expect more frequent occurrences of extremely poor ventilation conditions in winter over eastern China, which will reduce the efficacy of the emission control currently being implemented (34) but will provide a strong incentive for more stringent air pollutant and greenhouse gas emission reduction in China. Beijing won the bid to host the 2022 Winter Olympics Games. A repeat of the January 2013 weather condition will pose a much bigger challenge for regional emission control than the 2008 Summer Olympics Games because PM can be readily removed by rain events in summer and a faster summer ventilation rate also reduces the pollution contribution from regional accumulation compared to that in winter.


Observational data sets

The air pollution data for this study consist of four sources: (i) ground in situ PM10 monthly concentrations (2005 to 2015) retrieved from the daily air pollution index of five major cities [Zibo, Jining, Kaifeng, Pingdingshan, and Jinzhou (chosen because of the longest time series in these cities)] in the ECP region collected from the Ministry of Environmental Protection of the People’s Republic of China (, (ii) ground in situ PM2.5 concentrations in Beijing (2009 to 2015) collected from the Mission China air quality monitoring program of the U.S. Department of State (, (iii) the relative humidity–corrected meteorological Vil (1981 to 2013) (15) at the 45 plain ground sites (<300 m altitude) over the ECP region calculated from the Global Surface Summary of the Day (GSOD) database (version 8) provided by the National Climatic Data Center (NCDC), and (iv) the monthly AOD (2001 to 2015) at 550 nm derived from MODIS onboard the Aqua and Terra satellites (16). All these data sets, which are the currently available aerosol observations with the longest time series for the ECP region, have been widely used to study aerosol pollution in China (3538).

We examined two reanalysis data sets to evaluate ventilation conditions over the ECP region in January for the past three decades (fig. S5). The first one was the NCEP/NCAR reanalysis data (11) used in this study, and the second one, for cross-validation, was the ERA-Interim data (39) from the European Centre for Medium-Range Weather Forecasts. On the basis of the 1981 to 2015 reanalysis data, we computed ventilation indices with the following equationEmbedded Image(1)where Embedded Image is the normalized WSI (unitless) for the jth grid point of the ECP region in the ith year, Embedded Image is the monthly mean wind speed (in meters per second) at 1000 hPa for the jth grid point in the ith year derived from zonal and meridional winds of the reanalysis data, Embedded Image is the climatological monthly wind speed (in meters per second) for the jth grid point averaged from 1981 to 2010, and Embedded Image is the SD of wind speed (in meters per second) for the jth grid point from 1981 to 2010. Gridded atmospheric temperature gradient anomalies were normalized in the same manner based on the monthly potential temperature gradient between the fields at 925 and 1000 hPa. We then obtained WSI and ATGI for the ECP region by spatial averaging. It is necessary to have one ventilation index for the ECP region to simplify the interpretation of the multivariate statistical analysis results; thus, we calculated PPI as a synthetic meteorological index for each grid point to obtain the spatial distribution and then averaged it over the ECP region to obtain monthly time series using a weighted average of WSI and ATGIEmbedded Image(2)where r1 and r2 are the Pearson correlation coefficients of WSI (r1 = −0.73) and ATGI (r2 = 0.70), respectively, with in situ PM10 observations (table S2). The diagnostic ventilation indices (WSI, ATGI, and PPI) derived from the two reanalysis data sets agree well with each other, and the correlation between the two PPIs is 0.80 (fig. S5).

To investigate the association of climate factors to the ventilation condition, we collected multiple climate variables for the past three decades (1980 to 2015) in table S3. The first three meteorological indices were calculated based on the NCEP/NCAR reanalysis data (11) to describe the characteristics of the EAWM system (1820). The next two climate indices, Arctic Oscillation (internal atmospheric variability) (17) and ENSO (23), were collected from the Climate Prediction Center of the National Oceanic and Atmospheric Administration (NOAA). The last two cryospheric forcing factors were SIC (Fig. 1A) in the preceding autumn from the Met Office Hadley Centre (HadISST) (21) and boreal Eurasia SCE (Fig. 1A) in early winter from the Global Snow Lab at Rutgers University (22). The cryospheric indices, SIC and SCE, were normalized in the same manner as WSI and ATGI. We first averaged Arctic SICs within the Arctic Circle (66.6°N; Fig. 1A) in the preceding autumn and early winter seasons (August to November) and Eurasian SCE over the boreal region (60°E to 150°E, 40°N to 75°N; Fig. 1A) in early winter (October and November) for each year (1980 to 2014). We then normalized both variables with respect to their climatology (1981 to 2010) to obtain SIC and SCEEmbedded Image(3)where Xi is the ith year’s cryospheric variable such as Arctic SIC or Eurasian SCE, Xmean is the climatological average, Xstd is the SD for the same period, and Indexi is the normalized index in the ith year.

Statistical analysis

Statistical significance tests were used extensively throughout this study. We applied the moving block bootstrap method (12) to examine whether the January wind speed, the temperature gradient, and the Z850 daily data of 2013 are statistically different from the 30-year (1981 to 2010) climatological January data in Figs. 1 and 2. The moving block bootstrap method removes biases introduced by autocorrelation of the data of time length L or shorter (12). We first collected the daily NCEP/NCAR reanalysis data (11) and then regenerated the moving block bootstrap samples for each grid point with a block length of L = 5 days and a sampling size of n = 5000. The null hypothesis was that the 2013 data and the 30-year data were statistically from the same probability distribution with equal means. For those grid points with P values less than 0.01 (or 0.05), we rejected the null hypothesis and concluded that the values in 2013 over these areas were significantly different from the climatology at the 99% (or 95%) significance level. The same method was also applied to examine the significance of surface temperature anomalies in December 2012 in daily reanalysis data of fig. S4.

We used the standard bootstrapping method (12) to estimate the significance of correlation coefficients in table S2, PPI sensitivity responses in Fig. 3, and modeling surface air temperature responses in fig. S4 because all the time series used here were the monthly mean data of each year. For example, we estimated the correlation coefficients and their significance levels in table S2 on the basis of the empirical distribution of n = 5000 bootstraps.

We applied the MCA (12) to the PPI and Z850 fields to identify dominant circulation patterns affecting PPI over the ECP region. The MCA method performs a singular value decomposition of the covariance matrix of two variables to generate the coupled modes for the two variables separated in space and time dimensions. The temporal matrices are shown in Fig. 2A, and the spatial matrices are shown in Fig. 2 (B and C). The covariance of the two variables was maximized in the first mode. In our case, the first coupled MCA modes explained 33% of the covariance between the Z850 and PPI fields, 23% of the Z850 variance, and 35% of the PPI variance. These MCA modes and results remained consistent in a sensitivity test in which the 2013 data were excluded.

We examined the relationship between PPI and all the climate indices for the last 35 years listed in table S3 using the PCA. PCA is a dimension reduction method (12) to identify the major factors contributing to the variation of the variable of interest, which is PPI in this study. All the climate data (table S3) used in the PCA were first detrended. We added these data into a 35 × 7 matrix and computed the PCs. For attribution analysis, we applied the PC regression (PCR) method (40) to regress the detrended PPI against the PCs and examined their regression coefficientsEmbedded Image(4)where Zj(t) is the jth PC as a function of time, βj(j = 1, …, 7) is the corresponding regression coefficient, and Y(t) is the ECP regionally averaged PPI as a function of time.

Table S5 shows the regression coefficients and their P values for the F statistic of the hypotheses test to determine whether the corresponding coefficient is equal to zero or not. We followed the PC selection rule by Fekedulegn et al. (40) and found three major PCs (PC2, PC5, and PC6) that contributed significantly to PPI (P ≤ 0.05). Using the results from the PCR analysis, we found that these three PCs accounted for 53% of the PPI variance, whereas the inclusion of all PCs accounted for 57% of the variance.

The correlation coefficients of PCs with the detrended PPI in table S4 show similar results, that is, that PC2, PC5, and PC6 are the most significant PCs. The correlation coefficients of PCs with climate indices are also shown. Of particular interest to this study were the high correlations of SCE with PC5 (r = 0.83) and of SIC with PC6 (r = −0.80).

Note that the PCA was only applied to climate indices. The relationship between climate PC and detrended PPI was established using the PCR analysis. Therefore, the PCs explaining a large portion of the variance of climate indices were not necessarily correlated with PPI. For example, PC1 contributed most to the total variance of the matrix of climate indices, but it did not make significant contribution to the variance of detrended PPI. Interested readers are referred to a previous study (41) as an example of using the PCA to understand the effects of climate change on regional air quality.

We highlighted the PCA results in fig. S1. PC6, related to Arctic sea ice forcing, was the most important, contributing 29% to the 2013 PPI extreme over the ECP region, and explained 12% of the total variance in detrended PPI. PC5, related to Eurasian snow forcing and the Siberian High variability, was the second most important, contributing 13% to the 2013 PPI extreme, and explained 34% of the total variance. PC2, related to ENSO, was the third most important, contributing 8% to the 2013 PPI extreme, and explained 28% of the total variance. We added the PPI linear trend back to the contributions in fig. S1B. The PCA-reconstructed PPI explained most of the variations of the original PPI time series as well as the extreme in 2013, thereby placing the 2013 extreme in the context of the changes in the last 35 years. From 1981 to 2015, it was only in 2013 that the contributions of all three PCs are relatively large and, more importantly, all positive.

The two cryospheric forcing factors correlated strongly with PC5 and PC6, which had larger contributions to the 2013 extreme than the other PCs. To facilitate the subsequent analyses, we combined the two cryospheric forcing factors into a single normalized index, CFI, by weighted averaging SIC and SCE in a manner similar to the PPI formulationEmbedded Image(5)where r1 and r2 are the correlation coefficients of SIC (r1 = −0.43) and SCE (r2 = 0.64) with PPI, respectively.

Modeling experiments

The CESM (version 1.2.1) is a fully coupled global climate model maintained by the Climate and Global Dynamics Laboratory at the NCAR. We designed four experiments to investigate the impacts of cryospheric forcing on PPI over the ECP region using the component set of the CAM5 and the Community Land Model version 4.0 at a horizontal resolution of 1.9° × 2.5°. In the CTRL experiment, we conducted a 30-year simulation (with another 1-year simulation as spin-up) with prescribed climatological (averaged from 1981 to 2010) Arctic SIC (fig. S2) and SST from the Met Office Hadley Centre (HadISST) (21). In sensitivity simulations, we used January of each modeling year from the CTRL simulation as initial conditions for the corresponding sensitivity ensemble member and conducted a 14-month simulation with different cryospheric forcing. We analyzed the continuous December-January-February (DJF) data at the end of each sensitivity simulation to examine the seasonal impact of the cryospheric forcing. In the first sensitivity experiment (SENS1), we replaced the climatological SIC and SST data from August to November, which were used in the CTRL experiment, with observed cyclical (August 2012 to November 2012) HadISST data (21) over the Arctic region (fig. S2) following the method of Peings and Magnusdottir (42). In the second sensitivity experiment (SENS2), we introduced early winter snow forcing over boreal regions (60°E to 150°E, 40°N to 75°N; Fig. 1A and fig. S2) by perturbing October and November snowing rates of each ensemble simulation based on the observed SWE relative anomaly (24) in 2012 with the same climatological SIC/SST in the CTRL scenario. Therefore, we considered the snow forcing including the surface albedo effect and the insulation-related effect (for example, thermal conductivity and latent heat flux due to snow depth changes) in SENS2, both of which could cause significant local temperature response (43). The perturbation on modeling snow rates was performed using the following formulationEmbedded Image(6)where Snow0 is the default modeling snow rate without perturbations, Snow1 is the perturbed snow rate, and swefrac is the observation-based SWE fractional anomaly by comparing the 2012 anomalous SWE with the climatology (averaged from 1999 to 2010 because of limited data set availability; swefrac ≥ − 100 %; fig. S2). In the third sensitivity experiment (SENS3), we added both Arctic sea ice and boreal snow perturbations into the ensemble simulations to investigate their synergistic effects.

After numerical simulations, we calculated WSI, ATGI, and PPI over the ECP region for each scenario using the same method as that for reanalysis data. The CAM5 model uses a hybrid sigma-pressure vertical coordinate; thus, we extracted the modeling outputs at the nearest pressure levels to the reanalysis data in the calculation. The 30-year CTRL results were used as the climatological condition against which sensitivity results were compared. To evaluate the seasonal impact of cryospheric forcing, we had n = 90 ensemble members (DJF) for each sensitivity simulation in Fig. 3 and fig. S4.

To obtain the intensity of the dominant MCA modes in simulated Z850 and PPI fields, we projected the modeling fields onto the spatial patterns of the first MCA modes identified in the reanalysis data (Fig. 2, B and C) and estimated the MCA mode intensity of each ensemble member. The density distributions of the first pairwise MCA modes in fig. S3 were estimated by the two-dimensional kernel density estimation (kde2d) function included in the math package of R. By comparing density distributions of the sensitivity experiments to the control experiment, we evaluated the response sensitivities of both geopotential height and PPI fields to the specific cryospheric forcing.

Finally, we examined the ensemble results of the CESM1-CAM5 model (with three ensemble members) as well as its subset CCSM4 (with six ensemble members) in the CMIP5 project (13) based on their historical simulations (1980 to 2005) and future projections (2006 to 2100) under a medium mitigation scenario (RCP4.5) (table S6). We chose these two models because of their similarity with the model that we used in our sensitivity experiments. All the indices, including the cryospheric forcing in early seasons and PPI in January, were estimated in the same manner as described before. We then compared the correlation coefficients between the simulated cryospheric forcing and PPI with those based on reanalysis data to evaluate the modeling capability to reproduce the observed relationships. The differences between CMIP5 simulations and observations are discussed in Results.


Supplementary material for this article is available at

fig. S1. PCA decomposition and reconstruction of PPI.

fig. S2. Cryosphere forcing specifications used in the CESM numerical experiments.

fig. S3. Sensitivity response of the first coupled PPI and Z850 modes to cryospheric forcing.

fig. S4. Comparison of surface air temperature in December between reanalysis data and numerical experiments.

fig. S5. Time series of monthly WSI, ATGI, and PPI over the ECP region for January.

table S1. A list of acronyms used in this study.

table S2. Correlations among ventilation indices, PM observations, and cryospheric forcing factors.

table S3. Climate and synoptic weather indices in the PCA.

table S4. Correlation statistics among PCs and climate indices.

table S5. PCR coefficients of the detrended PPI onto PC.

table S6. The CMIP5 ensemble simulations used in this study.

Reference (44)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge high-performance computing support from Yellowstone ( provided by NCAR’s Computational and Information Systems Laboratory (CISL) and sponsored by the NSF. We thank the NOAA/Office of Oceanic and Atmospheric Research/Earth System Research Laboratory Physical Sciences Division (Boulder, CO, USA) for providing NCEP reanalysis data on the website at, the Met Office Hadley Centre for the HadISST data set, the NCDC for the GSOD data set, and the Global Snow Lab at Rutgers University for gridded SCE data. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in table S6 of this paper) for producing their model output and making it available. For CMIP, the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provided coordinating support and led the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We also thank Y. Peings, X. Li, and Z. Xie for helpful discussion on statistical data analysis and climate sensitivity experiments. Funding: This work was supported by the NSF Atmospheric Chemistry program and the U.S. Environmental Protection Agency (EPA) Science to Achieve Results programs through grant RD-83520401. It has not been subjected to any EPA review and therefore does not necessarily reflect the views of the Agency, and no official endorsement should be inferred. Author contributions: Y. Zou and Y.W. designed the analyses and modeling experiments. Y. Zou performed numerical simulations and prepared all figures. Y. Zou and Y.W. analyzed the data and wrote the paper. All authors discussed the data analysis and reviewed the manuscript text. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data are archived on the High Performance Storage Systems of the Yellowstone system managed by the CISL of NCAR and may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article