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.

Keywords

( 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.

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.

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.

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 PM 2.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 PM 2.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 ).

( 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).

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.

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.

( 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 ( ) 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 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 ( ). 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).

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).

( 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.

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.

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 ( 18 – 20 ), 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 ).

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 PM 10 over the ECP region, PM 2.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 PM 10 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 PM 10 (r = 0.92), PM 2.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 ).

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.

The winter cold extremes and heavy snowfalls in northern high latitudes are partially attributable to Arctic sea ice loss ( 28 – 30 ). 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.

MATERIALS AND METHODS

Observational data sets The air pollution data for this study consist of four sources: (i) ground in situ PM 10 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 (http://datacenter.mep.gov.cn/index), (ii) ground in situ PM 2.5 concentrations in Beijing (2009 to 2015) collected from the Mission China air quality monitoring program of the U.S. Department of State (www.stateair.net/web/post/1/1.html), (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 (35–38). 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 equation (1)where is the normalized WSI (unitless) for the jth grid point of the ECP region in the ith year, 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, is the climatological monthly wind speed (in meters per second) for the jth grid point averaged from 1981 to 2010, and 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 ATGI (2)where r 1 and r 2 are the Pearson correlation coefficients of WSI (r 1 = −0.73) and ATGI (r 2 = 0.70), respectively, with in situ PM 10 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 (18–20). 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 SCE (3)where X i is the ith year’s cryospheric variable such as Arctic SIC or Eurasian SCE, X mean is the climatological average, X std is the SD for the same period, and Index i 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 coefficients (4)where Z j (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 formulation (5)where r 1 and r 2 are the correlation coefficients of SIC (r 1 = −0.43) and SCE (r 2 = 0.64) with PPI, respectively.