A linear Markov model has been developed to predict sea ice concentration (SIC) in the pan-Arctic region at intraseasonal to seasonal time scales, which represents an original effort to use a reduced-dimension statistical model in forecasting Arctic sea ice year-round. The model was built to capture covariabilities in the atmosphere–ocean–sea ice system defined by SIC, sea surface temperature, and surface air temperature. Multivariate empirical orthogonal functions of these variables served as building blocks of the model. A series of model experiments were carried out to determine the model’s dimension. The predictive skill of the model was evaluated by anomaly correlation and root-mean-square errors in a cross-validated fashion . On average, the model is superior to the predictions by anomaly persistence, damped anomaly persistence, and climatology. The model shows good skill in predicting SIC anomalies within the Arctic basin during summer and fall. Long-term trends partially contribute to the model skill. However, the model still beats the anomaly persistence for all targeted seasons after linear trends are removed. In winter and spring, the predictability is found only in the seasonal ice zone. The model has higher anomaly correlation in the Atlantic sector than in the Pacific sector. The model predicts well the interannual variability of sea ice extent (SIE) but underestimates its accelerated long-term decline, resulting in a systematic model bias. This model bias can be reduced by the constant or linear regression bias corrections, leading to an improved correlation skill of 0.92 by the regression bias correction for the 2-month-lead September SIE prediction.
The rapid summer Arctic sea ice retreat has not only been an icon of climate change but has also created more commercial opportunities in the newly opened Arctic waters, such as shipping and oil drilling (Eicken 2013). However, lower summer sea ice cover comes with larger ice variability (Goosse et al. 2009), causing tremendous difficulties in planning, and even threatening, commercial operations in the Arctic. Therefore, skillful sea ice prediction will become a needed service in the warming Arctic environment. The scientific questions that we are currently facing are the following: 1) To what extent are sea ice anomalies from year to year predictable, and 2) what skill do stochastic models have beyond the information in the long-term trend, anomaly persistence, or climatology?
Sea ice is a fully interactive component of the climate system. Sea ice concentration (SIC), extent, and thickness can vary on many different time scales, from intraseasonal, seasonal, interannual, to decadal and longer. Many studies have shown that sea ice variations are predominantly caused by atmospheric circulation changes (Fang and Wallace 1994; Deser et al. 2000; Pfirman et al. 2004; Wu and Zhang 2010). Proshutinsky and Johnson (1997) suggest that the Arctic winter sea ice variability is forced by changes in the location and intensity of the Icelandic low and the Siberian high. Francis and Hunter (2007) reveal that the Bering Sea ice variability is influenced mainly by anomalies in easterly winds associated with the Aleutian low, whereas the Barents Sea ice variability is driven primarily by sea surface temperature (SST) and local meridional winds. Near the Alaska coastal region, the rapid sea ice decay is related to downward solar radiation and mechanical breakup due to winds and ocean currents following the thermal weakening of sea ice, which leads to sea ice predictability up to two weeks in advance (Petrich et al. 2012). A number of studies found that Arctic Oscillation variability affects Arctic sea ice on a broad range of time scales, ranging from seasonal to decadal (Deser et al. 2000; Rigor et al. 2002; Serreze et al. 2003; Rigor and Wallace 2004). Other low-frequency large-scale climate patterns (e.g., NAO, PNA, ENSO, Arctic dipole) also exert a strong impact on Arctic sea ice (Mysak et al. 1996; Mysak and Venegas 1998; Moritz et al. 2002; L’Heureux et al. 2008; Wang et al. 2009). These findings suggest that sea ice predictability likely lies in the covariability of the polar climate system. They imply that Arctic sea ice could be better predicted by empirical models when developed in the multivariate space.
Sea ice itself is a relatively slow-varying component in the atmosphere–ocean–sea ice coupled system, and SIC anomalies usually persist for a few seasons (Chapman and Walsh 1991; Yuan 2004). These characteristics suggest that sea ice has a certain level of natural predictability. A good example is the seasonal forecast of the Antarctic sea ice by a linear Markov model (Chen and Yuan 2004). The success of the model is attributed to the dominance of the Antarctic climate variability by a few distinctive modes in the atmosphere–ocean–sea ice coupled system, and to the model’s ability to pick up these modes. In the Arctic, a number of studies suggest that the summer sea ice may be predictable at a seasonal time scale because of links between summer sea ice condition and early season climate indices (Walsh and Zwally 1990; McPhee et al. 1998; Maslanik et al. 1999; Kwok 2000; Zhang et al. 2000; Drobot and Maslanik 2002). Chapman and Walsh (1991) found that an Arctic sea ice anomaly can persist for several months but a simple analog model has a meaningful skill of only one month. Drobot and Maslanik (2002) were able to predict sea ice anomaly in the Beaufort Sea a few months in advance by using a linear regression model. A few similar regression forecast models have been developed in the Beaufort and Chukchi Seas (Drobot 2003), in the Hudson Strait (Gough and Houser 2005), and for the pan-Arctic summer sea ice minimum extent (Drobot et al. 2006). This type of sea ice prediction model commonly uses the multiyear ice gradient, ice concentration, air temperature, and climate indices like the NAO index as predictors. Gough and Houser (2005) also found that the oceanic heat gain or loss during previous seasons, which is determined by the length of ice-free or ice-covered season, is fundamentally important to the regional forecast.
While the majority of statistical models were developed mostly for isolated regions or specific seasons, coupled climate models, on the other hand, are able to predict sea ice concentration and ice thickness at pan-Arctic grid points and for all seasons. Many dynamic models show some predictive skill at seasonal to interannual time scales (Collins 2002; Koenigk and Mikolajewicz 2009; Blanchard-Wrigglesworth et al. 2011a,b; Chevallier et al. 2013; Guemas et al. 2014; Tietsche et al. 2014). For predicting September Arctic sea ice extent (SIE), Lindsay et al. (2008) found that sea ice concentration is the best predictor in the first two months, while the ocean temperature at 200–270 m is more important for longer lead-time prediction in their ice–ocean coupled model. In addition, studies revealed that sea ice anomalies in melting seasons reemerge in freezing seasons in multiple models, which is caused by the persistence of SST anomalies near the ice edge (Blanchard-Wrigglesworth et al. 2011a; Day et al. 2014; Bushuk et al. 2014, 2015). Chevallier et al. (2013) found that the Centre National de Recherches Météorologiques Coupled Global Climate Model has significant skill in predicting the Arctic SIE in September even after the linear trend was removed. On the other hand, the prediction skill for the Arctic sea ice extent was often masked by the effect of long-term trends in other coupled climate models (Lindsay et al. 2008; Sigmond et al. 2013; Msadek et al. 2014). Moreover, regional sea ice prediction skill from dynamic models was evaluated far less thoroughly at intraseasonal and seasonal time scales.
Although there are limitations to sea ice predictability owing to inherent chaotic atmospheric variability in the Arctic climate system (Holland et al. 2010; Blanchard-Wrigglesworth et al. 2011b), the above studies suggest that it is feasible to predict Arctic sea ice at intraseasonal to seasonal time scales. Despite great efforts toward predicting the September total SIE through the sea ice prediction network (http://www.arcus.org/sipn), seasonal or intraseasonal predictions for either the total ice extent or ice concentration at grid points are very limited in other months in the pan-Arctic region. Skillful sea ice forecast in the Arctic is still lacking even though the demand for such prediction has grown tremendously in recent decades. Because of its cost-effectiveness, we are motivated to develop more advanced statistical models for sea ice prediction in the pan-Arctic region.
In this study, we construct a linear Markov model for monthly SIC at pan-Arctic grid points for all seasons. The sea ice predictability is assessed at different locations and over all seasons. Through development of the model, we isolate the relative importance of different processes that dictate sea ice variability. In addition, we generate SIE forecasts based on predicted SIC and evaluate SIE prediction skill, with emphasis on the September SIE.
2. Data and methodology
In this study, we primarily use SIC data from the GSFC bootstrap SMMR–SSM/I version-2 quasi-daily time series from 1979 to 2013 (Comiso 2000). This 35-yr time series can be obtained from the EOS Distributed Active Archive Center (DAAC) at the National Snow and Ice Data Center (NSIDC; http://nsidc.org). Monthly sea ice concentrations from original 25 km × 25 km grids were bin-averaged into 2° × 0.5° longitude–latitude grids. Atmospheric variables are taken from NCEP–DOE AMIP-II reanalysis data with spatial resolution of 2.5° × 2.5°. This dataset has been improved from the original NCEP–NCAR reanalysis (Kalnay et al. 1996; Kistler et al. 2001) by fixing errors and updating the parameterization of the physical process in the model (Kanamitsu et al. 2002; http://www.cpc.ncep.noaa.gov/products/wesley/reanalysis2/kana/reanl2-1.htm). In addition, we use SST from Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) monthly data. The SST data were reconstructed from in situ observations and satellite measurements from various sources with a spatial resolution of 1° × 1° longitude–latitude (Rayner et al. 2003).
To isolate which oceanic and atmospheric variables are closely related to sea ice variability, we calculated correlations between SIC and the following: SST, surface air temperature (SAT), pressure and winds at sea level, and geopotential heights and winds at 500-hPa and 300-hPa levels, at collocated grid points. On average, SIC is highly related to SST and SAT, followed by geopotential height at 300 hPa, and then geopotential height at 500 hPa, as well as winds at these levels. The correlations with sea level pressure and surface winds are relatively small. Because of the barotropic nature of the polar troposphere (Ting 1994; Chen 2005), only one pressure level in the troposphere is chosen. Therefore, SIC, SST, SAT, and geopotential height and winds at 300 hPa are chosen to form the initial multivariate space for capturing the predictable variability in the atmosphere–ocean–sea ice system. A linear Markov model is then constructed in the multivariate EOF (MEOF) space. The base functions of the model’s spatial dependence consist of the eigenvectors from the MEOF of these six variables, while the temporal evolution of the model is a Markov process with its transition functions determined from the corresponding principal components (PCs). We use only several leading MEOF modes, which greatly reduce model space and filter out unpredictable small-scale features. Similar Markov models have been used to predict Antarctic sea ice (Chen and Yuan 2004), El Niño–Southern Oscillation (ENSO; Blumenthal 1991; Xue et al. 1994; Canizares et al. 2001; Xue et al. 2000), and East Asian monsoon (Wu et al. 2013). Here, we follow Chen and Yuan (2004) and construct a similar model for seasonal prediction of Arctic sea ice concentrations.
First, we create anomaly time series of SIC, SST, and atmospheric variables from 1979 to 2012 by subtracting climatologies of the same period from monthly mean data and normalizing all anomaly series. The normalization is applied to time series at grid points but does not account for the differences in the numbers of grid points for each variable. The number of grid points is 6728 for SIC, 5040 for SST, and 3024 for each atmospheric variable. To emphasize sea ice variability in the model construction, we weight SIC, SST, SAT, 300-hPa height, and u and υ at the 300-hPa level by 2, 1, 1, 1, 1, and 1, respectively, although the final model skill is not very sensitive to the weight. The weighted variables are stacked up into a single matrix (n, m), where n is the number of grid points of all fields and m is the length of time series. We then decomposed into eigenvectors (spatial pattern) and their corresponding PCs (time series) :
where the columns of are orthogonal and the columns of are orthogonal; the superscript T denotes matrix transposition. By truncating (1) to the first several modes, we greatly reduce the model space. The Markov model is computed using the single-step correlation matrix—that is, a transition matrix that satisfies the following linear relationship:
where i denotes the ith month and ei is the error in the model fit. Multiplying (2) with and averaging over time yields the following:
where angle brackets indicate time averaging. There should be no correlation between and for the best model fit, so can be calculated as follows:
Since sea ice behaves dramatically differently in different seasons, needs to be considered season dependent; thus (4) is actually applied to 12 subsets of PCs to obtain different transition matrices for each of the 12 calendar months. When the model is constructed, the sea ice prediction will be made through following four steps: 1) the PCs corresponding to the initial anomalous climate conditions are calculated by projecting observations to the eigenvectors , 2) the predictions of the PCs are made at increasing lead times by successively applying the transition matrices, 3) the predicted PCs are combined with the respective eigenvectors to forecast sea ice concentration anomalies, and 4) the predicted sea ice anomaly field is finally added into the climatology to produce the final prediction.
To determine model variables that need to be used in the model, we calculate the model’s hindcast skill measured by anomaly correlation (AC) and root-mean-square errors (RMSEs) between predictions and observations. For hindcast experiments, the model is initialized from each month of the 34 yr of time series, and 12-month predictions are made from each initialized month. The hindcast skill in this study is equivalent to a reconstructed time series skill (or in-sample skill). The model variables are then determined by trial and error based on hindcast skill averaged over 34 years. The MEOF is recalculated each time when a new variable is added. When we determine the model dimension, we evaluate the model skill in a cross-validation fashion. To avoid artificial skill, we subtract one year from PCs and recalculate the transition matrix in (4); then we make 12-month predictions for that year. Here we subtract one year of data from PCs instead of the original time series because the spatial patterns of leading modes remain essentially the same when one year of data is taken out from 34 yr of time series. The process is repeated for each year of the time series. In this process, the MEOF is not recalculated and the time series of PCs are reduced to 33 years. Therefore, with slightly shorter time series, the model prediction skill is assessed over 34 one-year predictions. The cross-validated evaluation is stricter than the hindcast assessment since the predicted year information was not included in the Markov model construction.
In addition, the hindcast skill and cross-validated skill are also compared to the skill of the SIC anomaly persistence, damped anomaly persistence, and climatology predictions. To be consistent with the hindcast skill, the persistence function is constructed as follows: the autocorrelations of SIC anomalies at each grid point are computed from 1- to 12-month lag starting from each month of the 34-yr time series and then are averaged as a function of lag month. Damped anomaly persistence provides a forecast that amplitudes of the anomalies are assumed to reduce in time exponentially at a time scale of the local autocorrelation. It is defined as , where is the initial anomaly, t is the lead time, and the slope α is determined by the system’s autocorrelation (Griffies and Bryan 1997). The damped persistence is also calculated and averaged in the same fashion as the anomaly persistence.
3. Arctic sea ice variability and MEOF analysis
Sea ice concentration and its total extent exhibit significant declining trends in all seasons (Parkinson and Cavalieri 2012; Simmonds 2015). During the period from 1979 to 2012, ice concentration trends vary greatly in space and are particularly large in the Barents Sea, in the Kara Sea, and north of the Chukchi Sea (Fig. 1a). The variability of ice concentration measured by standard deviation of anomaly is also high in these areas with large linear trends (Fig. 1b). Because the trends are significant parts of the total variability, we decided to retain the trends in all the time series. The MEOF analyses were carried out on the normalized, weighted, and combined data matrix. Figure 2 shows the eigenvectors of the first three MEOF modes of all model variables, which account for about 11%, 6%, and 4% of the total variance, respectively. Their corresponding PCs are displayed in Fig. 3. The first mode clearly represents the long-term trends of SIC declining and SST and SAT warming. The first-mode SIC spatial pattern (Fig. 2) resembles the sea ice long-term trends in Fig. 1a. The trend is also clearly evidenced in PC1 (Fig. 3a). These SIC and temperature trends accompany the increasing of geopotential heights in the mid–high latitudes of the Northern Hemisphere, particularly in the Siberian high and the subtropical high in eastern North America (Fig. 2). The second mode represents a typical quadruple pattern of sea ice anomalies in winter (Yang and Yuan 2014); the same phase of ice anomalies is found in the Barents and Kara Seas and Sea of Okhotsk, while the opposite phase of ice anomalies is found between the Barents and Kara Seas and Baffin Bay and between the Okhotsk Sea and Bering Sea (Fig. 2e). Surface temperature anomalies are consistent with the sea ice anomalies, while geopotential height anomalies resemble a positive AO pattern (Fig. 2h). This mode also contains summer sea ice variability in the Beaufort, Chukchi, East Siberian, and northern Barents Seas. The SIC summer variability is weakly anticorrelated to the winter variability as revealed by −0.2 autocorrelation of PC2 at 6-month lag. The physical process causing the relationship is yet to be understood. The interannual variability dominates the mode-2 temporal variability. It is noted that the amplitude of temporal variability is rather low between 1995 and 2005 for these three modes. The reason for the low-variability period remains unknown, except the trend is small in the midpart of time series for PC1. The eigenvalues from this MEOF analysis (see Fig. S1a in supplementary material file JCLI-D-15-0858s1) reveal that the first four modes are significantly different from their neighboring modes based on the North test (North et al. 1982). Clearly, the three modes presented in Figs. 2 and 3 represent the sea ice–related large-scale climate variability in the Northern Hemisphere.
4. Model experiments
The model experiments were first conducted to investigate contributions from each variable and to test model sensitivity to the number of MEOF modes included in the model. After selecting the variables to include and determining the number of MEOF modes for the model, we then evaluated the model skill in a cross-validated fashion. Finally, the SIE forecasts were made and the skill of SIE predictions was assessed with and without linear trends.
a. Hindcast experiments for determining model variables
Here we examine whether the six variables in the initial multivariate space are necessary for achieving the best prediction skill. The MEOF analysis was first applied to SIC alone. Then SAT, SST, 300-hPa geopotential height, and winds at the 300-hPa level were combined with SIC separately. Finally, all six variables were added into the data matrix. The MEOF analysis was repeated to each constructed data matrix. Using the PCs from each set of MEOF analyses, the Markov models were developed using (2)–(4), and hindcast predictions were carried out for each month of the time series with 1- to 12-month lead times. Figures 4 and 5 show the hindcast skill measured by AC and RMSE, respectively, in the six areas identified in Fig. 1b. Both AC and RMSE were averaged at grid points over all initial months before the areal mean was calculated. The averaged skill is not overly sensitive to the size of averaging area in each sea. First, the model skill significantly exceeds the SIC anomaly persistence in terms of AC and drastically reduces RMSE from the anomaly persistence, regardless of what variables are included. Second, prediction skill is highly variable throughout the Arctic. The skill is quite high in the Chukchi Sea, above 0.6 even at a 12-month lead. It is followed by the skill in the East Siberian Sea, Barents Sea and Baffin Bay. The skill in the Sea of Okhotsk and Bering Sea is relatively low. Third, different regions need different combinations of model variables. For AC calculation with a sample size of 396 for any given lead time, the 0.1 difference between two ACs (both >0.5) would be significant at the 95% confidence level, based on a two-tailed Fisher r-to-z transformation. The model built on the data matrix of SIC, SAT, and SST performs better in most areas except in the Chukchi Sea and Bering Sea, where geopotential height and winds contribute most to the skill in 1- to 4-month-lead predictions. It reflects that dynamic processes are more important than the thermodynamic process to the short-term SIC predictability in these regions. On the other hand, the thermodynamic process dominates the model skill in lead times longer than four months in the Chukchi Sea. The RMSE show similar characteristics revealed by AC. It is worth mentioning that the skill in the Arctic peripheral seas mainly reflects the model performance in winter and spring, while the skill in marginal seas within the Arctic basin reflects the skill in summer and fall in Figs. 4 and 5. The seasonally averaged skill is quite close to the annually averaged skill since zero anomaly data points (such as summer SIC in the Bering Sea and winter SIC in the Chukchi Sea) do not contribute much to the correlation.
Since surface temperatures seem to contribute most to the model skill, we ran another set of experiments to isolate the contributions from SST and SAT. SST clearly plays a more important role than SAT (not shown). Including both SST and SAT slightly outperforms the case without SAT; thus, SAT is a minor contributor but can still contribute to the skill. Since most areas and lead times we tested show that the combination of SIC, SST, and SAT produces better prediction skill (Figs. 4 and 5), we decide to use these three variables to construct the model.
Although sea ice concentration is observable from space, it alone does not reflect all sea ice dynamic behavior. We also conducted model experiments that include sea ice thickness (SIT). The ice thickness data were taken from a pan-Arctic ice ocean modeling and assimilation system (PIOMAS; Zhang and Rothrock 2003). The Markov model quickly gains skill when it includes SIT in the Barents, East Siberian, and Chukchi Seas (Figs. 6b–d). Including atmospheric variables and SST only slightly enhances skill there. In the Baffin Bay, ice thickness dominates the prediction skill for all lead months. Adding atmosphere variables and SST only reduces the skill (Fig. 6a). Surprisingly, adding sea ice thickness to the model has a negative contribution to the prediction skill in the Bering Sea and the Sea of Okhotsk (Figs. 6e,f). Why ice thickness reduces prediction skill in these regions remains unknown. One possible reason is that PIOMAS ice thickness has larger errors in the Bering Sea and the Sea of Okhotsk. The Markov model skill is likely sensitive to the quality of ice thickness data. This hypothesis is difficult to verify owing to the lack of observations to evaluate ice thickness in PIOMAS and other oceanic reanalysis products (Chevallier et al. 2016), particularly in the Bering Sea and the Sea of Okhotsk. On the other hand, by using the oceanic and atmospheric variables presented here, the model captures the majority of predictable variance without having access to the ice thickness information. It reflects the effectiveness of building a Markov model in the MEOF space. Therefore, we do not rely on ice thickness derived from PIOMAS in our model development.
It is worth mentioning two SIC characteristics presented in the persistence prediction. First, both persistence prediction skill and persistence RMSE in Figs. 4 and 5 show low predictability and large errors at around a 6-month lead time. It reflects the fact that summer SIC anomalies in the Arctic basin cannot predict winter ice condition by the persistence model because of low SIC variability due to the total freeze-up in winter. Similarly, winter SIC anomalies in the Arctic peripheral seas will not lead to a correct prediction in summer by the persistence model owing to complete sea ice melt in summer. In such a case, the persistence and damped persistence models are worse than climatology predictions (Fig. 5). On the other hand, the Markov model does not have such shortcomings at the 6-month-lead prediction. Second, the skill of the SIC persistence in some locations increases after reaching its minimum (see Figs. 4 and 6), reflecting the reemergence of SIC anomalies suggested in early studies (Blanchard-Wrigglesworth et al. 2011a; Day et al. 2014; Bushuk et al. 2014, 2015). The reemergence of SIC anomalies occurs at 6-month (the SIC anomaly in the melting season reappears during the freezing season) and year-to-year time scales. Since the SIC persistence presented here is averaged from all seasons, the 6-month reemergence is smoothed out but the annual reemergence stands out clearly. This reemergence characteristic partially contributes to the long-lead predictability of SIC anomalies.
b. Cross-validated experiments for determining model truncation
Eigenvalues of the MEOF analysis with SIC, SST, and SAT also suggest that only the first four MEOF modes are well separated from their neighboring modes (Fig. S1b in supplementary material file JCLI-D-15-0858s1). The question is whether these four modes capture most predictable variability in the pan-Arctic region. To determine a reasonable model truncation, we run cross-validated experiments with the model retaining 1 to 13 modes, respectively. The experiments are repeated over 34 years for a robust assessment as described in the section 2. The demand on the number of modes in the Markov model varies greatly across the Arctic (Figs. S2 and S3), so we decide to choose a configuration that better fits most areas and seasons. Therefore, the model skill is evaluated using the following three globally averaged measures. First, we calculated the mean AC for all grid points that have significant correlations. Second, we calculated the percentage of grid points that have significant correlations in the latitude band of 50°–80°N for winter and spring, and 75°–90°N for summer and fall. The purpose is not to compare the differences between winter and summer but to identify skill differences induced by retaining different number of modes within these seasons. Third, we compute the mean RMSE in the areas where ACs are significant. These evaluations focused on the areas that exhibit significant sea ice predictability. We consider that the model is overfitted when the AC remains the same (or reduced) but RMSE is increased. The model with 1 to 5 modes clearly has lower skill than that with 7 to 13 modes in all seasons (Figs. 7a,b). In summer and fall, the model with one to five modes also has a relatively low percentage of significant correlation and higher RMSE (Figs. 7c,e). To our surprise, the model with one to five modes has a higher percentage of grid points with significant correlation and lower RMSE in winter and spring, although skill is lower (Figs. 7b,d,f). On the other hand, the mean skill of the model experiments with 7 to 13 modes is quite similar, in which the model with 11 and 13 modes has slightly better skill and higher percentage of significant correlations up to a 4-month lead time as well as lower RMSE in summer and fall. However, the model with 13 modes shows a clear lower percentage of significant correlation and higher RMSE in winter and spring (Figs. 7d,f). It indicates that including 13 modes does not bring more predictability but likely introduces more unpredictable noise. To avoid overfitting the model, we decide to include 11 modes in the model. Although MEOF modes 5 to 11 are not significantly separated from their neighboring modes, they represent some regional characteristics of ice behaviors in different seasons. The prediction skill of the Arctic SIE will further confirm that including a high number of MEOF modes would capture more natural variability in sea ice (see next section).
Figure 7 also reveals characteristics of the pan-Arctic sea ice predictability in this model. On average, the prediction skill is relatively high in summer and fall, ranging from 0.7 at a 1-month lead to 0.4 at a 12-month lead. The prediction skill for winter and spring is about 0.1 ~ 0.2 lower than that in summer and fall. On the other hand, the high prediction skill comes with relatively high RMSE in summer and fall (cf. Fig. 7e and Fig. 7f). The predictability of Arctic sea ice clearly exhibits seasonal and regional characteristics, which need to be explored further.
c. Assessment of model skill
Once the final model configuration is decided (the model with three variables and 11 MEOF modes), the MEOF of the three variables remains fixed. The model prediction skill was evaluated at all grid points and for all seasons in the cross-validated fashion for the 34-yr period. For winter (DJF) prediction, high skill is concentrated in the Atlantic sector of the Arctic marginal seas and thier peripheral seas: the Barents Sea, Kara Sea, Baffin Bay, and Labrador Sea. The skill is only slightly reduced from a 3-month to 9-month lead. The Sea of Okhotsk and the Bering Sea have relatively low prediction skill. This is likely because dominant coupled climate variability occurs in the Atlantic sector and covariability with smaller amplitudes unique to the Pacific sector is not picked up by the leading modes. Moreover, the geopotential height and wind information that contribute to the model skill in the Bering Sea were not included in the model. The skill in the central Arctic basin is also low because the area is nearly 100% covered by sea ice in winter (Fig. 8). Prediction skill in spring has a similar pattern to that in winter but with lower predictability. In summer (JJA), sea ice in the Arctic peripheral seas is melted. All prediction skill is concentrated within the Arctic basin: the 3-month-lead prediction has skill above 0.6 for most of the central basin. Reduced skill is found near the marginal seas, particularly for the longer lead time. The prediction skill in fall has a similar spatial pattern as that in summer but with higher skill, particularly for 9-month-lead predictions. This could be related to the Arctic sea ice anomaly reemergence from spring to fall due to the oceanic memory (Blanchard-Wrigglesworth et al. 2011a; Bushuk et al. 2014, 2015). In general, the model has higher prediction skill for summer and fall than for winter and spring, while skill is the lowest in spring.
In contrast to the higher ACs in the Atlantic sector than the Pacific sector, the RMSEs are more uniformly distributed around the Arctic, slightly higher in the Barents Sea and Bering Sea in winter and spring (cf. Fig. 8 to Fig. 9). The magnitudes of the RMSE stay at the same level from a 3-month to 9-month lead. The pattern and magnitudes of RMSE for summer and fall predictions are also quite similar, although the marginal seas have relatively larger errors than the central basin (Fig. 9). Occasionally, RMSE can reach above 25% of the ice concentration, which is a quite substantial error. However, RMSE is consistently smaller than the standard deviation of ice concentration anomalies across the Arctic (Fig. S4 in supplementary material file JCLI-D-15-0858s1).
In addition, the model performance is further evaluated against three metrics: anomaly persistence, damped anomaly persistence, and climatology. The damped anomaly persistence and anomaly persistence behave similarly in terms of AC but produce different RMSE (Wang et al. 2016). So we evaluate the model prediction against these metrics through RMSE. Averaged over pan-Arctic grid points and over 34-yr model experiments, the mean RMSE of Markov model predictions of SIC does not change much with lead time but is consistently smaller than that produced by anomaly persistence, damped anomaly persistence, and climatology predictions for all four seasons (Fig. 10). The model’s mean RMSE is larger in summer and fall than in winter and spring, reflecting larger SIC variability in summer and fall. Even though the mean RMSEs from Markov model predictions are only slightly smaller than climatology predictions in spring (Fig. 10), the differences between these two sets of RMSEs are statistically significant owing to the large sample size (357 390).
To assess the performance of the Markov model relative to dynamic forecast systems, we calculated the skill from the NCEP Climate Forecast System, version 2 (CFSv2; Wang et al. 2013), which has similar performance compared to other dynamic forecast systems (Merryfield et al. 2013; Blanchard-Wrigglesworth et al. 2015). The Markov model appears superior in most areas that we examined in terms of prediction skill (Fig. 11). The skill of these two models is comparable only in the Sea of Okhotsk. It is worth mentioning that the hindcast skill in this statistical model is not calculated in the same way as the hindcast skill in the dynamic model. However, they are the available metrics for comparisons between different types of models. The model also shows an advantage over earlier statistical models, such as the analog model presented in Chapman and Walsh (1991), which cannot outperform persistence beyond a 1-month lead.
d. Prediction of Arctic sea ice extent
The Arctic SIE forecast in all months is calculated by summing all areas that have at least 15% sea ice concentration from the SIC predictions produced by cross-validated experiments and from observations. Here we use September SIE forecast to illustrate the model’s skill in predicting Arctic SIE. Figure 12 shows the 2-, 3-, and 4-month-lead September SIE predictions during 1980 to 2012 and the forward forecasts in 2013 and 2014. The model generally predicts the phase of SIE variability well but with weaker amplitudes of variability. It was able to predict the historical low in 2007 and 2012 and also captured the bounce back in 2013 and 2014 in the forward forecast (Figs. 12a–c). However, it did not capture the accelerated SIE decline, particularly in the last decade. The model skill of predicting September SIE measured by ACs and RMSE is listed in Table 1. These model assessments reflect its ability to predict anomaly phases (high AC) and its deficiency in accurately predicting magnitudes of anomalies. The model systematically overestimates the September SIE. One possible reason for this model behavior is that the model is developed in the multivariate space. Sea ice likely has the strongest trend among all oceanic and atmospheric variables. The trends in PCs of MEOF modes may not present the true strength of sea ice trends. To overcome this model deficiency, we consider two types of model bias corrections for SIE predictions.
The first-order or mean bias correction method can be achieved by subtracting a mean prediction error during 1980–2012 from all predictions. The bias correction reduces overestimation of SIE in the recent decade but produces underestimation of SIE in early years. This method does not change the model correlation skill but reduces approximately 30% of RMSE for all predictions at 2- to 4-month lead times. (Table 1).
The second approach of bias correction is to predict model errors. The forecast errors display a relatively good linear relationship with SIE anomalies in the initial months for 2- and 3-month-lead forecasts, although such a relationship vanishes in the 4-month-lead forecast (Figs. 12d–f). This indicates that SIE has limited duration of memory. The linear relationship allows us to conduct a bias correction using linear regression for the short-term predictions:
where SIEb represents the n-month-lead forecast bias (error) at the tth month and is expressed as the linear regression of the SIE anomaly SIEa at the (t − n)th month. The slope of the regression α represents the linear dependence on SIEa and the intercept β, the constant bias at zero anomaly. Such a bias correction does not use information beyond the lead time and is thus a forecast method for the SIE residual that has not been resolved by the linear Markov model for SIC. This method is most effective for the 2-month-lead forecast, which shows both enhanced correlations and reduced RMSE (Fig. 12 and Table 1). Specifically, the method reduces 62% of RMSE from the prediction without bias correction and removes an additional 41% of RMSE from the prediction with the constant bias correction. The skill for the 3-month-lead prediction was improved slightly from the results of prediction by constant bias correction. However, it does not show much improvement from the prediction with constant bias corrections for 4-month-lead predictions (Table 1). Because of a very low value of α (0.02) for the 4-month-lead prediction, the linear regression bias correction is equivalent to the constant bias correction. Both bias corrections improve model predictions mainly through reducing the model’s RMSE.
This linear regression bias correction is different than the methods presented in Kharin et al. (2012) and Fučkar et al. (2014). In Kharin et al. (2012), both predictions and observations can be expressed separately by a linear trend plus anomalies relative to the linear trend. The linear regression bias correction is then achieved by replacing the linear trend in predictions with the linear trend in observations. In Fučkar et al. (2014), predictions and observations were expressed separately by linear regressions of observed initial conditions plus anomalies with respect to the linear regressions. The regression bias correction is then achieved by replacing the linear regression in the model forecast with the linear regression in observations. For SIE prediction, these three types of linear regression bias correction produce similar results (Wang et al. 2016).
e. Contribution of linear trends to the prediction skill
In general, there are two types of predictability in a climate system: one is forced or derived from changing boundary conditions; another comes from initial conditions or internal variability (Lorenz 1975). In the Arctic, we consider that the long-term trend is forced by the climate change. Because of the rapid decline of Arctic sea ice, linear trends in sea ice account for a large portion of prediction skill in dynamic models (Lindsay et al. 2008; Blanchard-Wrigglesworth 2011b; Sigmond et al. 2013; Msadek et al. 2014). To isolate the contribution of sea ice trends to the prediction skill in the linear Markov model, we conducted postprediction analysis on the predictions from cross-validated experiments. Monthly trends at each grid point were removed from the predictions and observations. The skill of detrended prediction is evaluated by ACs of the residual series. On average, the trend removal reduces the skill by 0.16 (0.25) of AC for the 3-month (6 month) lead predictions and leaves 9-month predictions without any skill for all seasons (Fig. 13). Large skill reduction occurs north of the Chukchi Sea as well as in the East Siberian Sea and the Barents–Kara Seas, where SIC has large declining trends (Fig. 1a). However, the model still retains the predictive skill (AC > 0.6) in the central Arctic for 3-month-lead summer and fall predictions and in the Canadian basin for 6-month-lead fall predictions (Fig. 13). In winter and spring, the peripheral seas in the Atlantic sector also lose predictive skill after the trend removal (cf. Fig. 8 and Fig. 13) but retain predictive skill for 2-month-lead winter prediction (not shown). Although the trend removal largely reduces SIC prediction skill, the residual skill still beats the detrended persistence prediction in most areas in the Arctic basin for summer (fall) predictions up to 6-month (9 month) lead time, as well as in Arctic peripheral seas for 3-month-lead winter and spring predictions (Fig. S5 in supplementary material file JCLI-D-15-0858s1). This means that the model captures certain SIC internal variability.
Furthermore, we examine the trend’s contribution to the prediction skill of Arctic SIE. The model has skillful predictions of Arctic SIE from November to January for most lead times and high skill predictions from July to September up to 8-month lead times (Fig. 14a). Then we remove the monthly trends from SIE predictions and observations. As shown in section 4d, the Markov model underestimates the strength of the linear trend in September SIE. The model behaves similarly in other months. Removal of monthly linear trends independently from the predictions and observations actually improves the prediction skill for all seasons, which shows the ability of the model in predicting year-to-year variability of Arctic SIE. After the trend removal, the model is skillful (AC > 0.6) in predicting SIE from September to February but has low predictability from March to June (Figs. 14a,b). The improvements are especially large in June and October because the model produces larger errors in predicting the trends in these two months. Before the trend removal, the model’s SIE prediction can only beat anomaly persistence in isolated months with limited lead times (Fig. 14c). After the trend removal, the model can beat the persistence in most months and most lead times (Fig. 14d). The SIE forecast by the Markov model has smaller RMSEs than that by climatology prediction, except for spring forecasts with trends or without trends (Fig. S7 in supplementary material file JCLI-D-15-0858s1). Compared to the Canadian seasonal and interannual prediction system (Fig. 2 in Sigmond et al. 2013), the Markov model is more capable of capturing natural variability during all seasons but is less capable of correctly predicting long-term trends in SIE.
On the other hand, if the Markov model only retains four independent MEOF modes (Fig. S1b), the model loses overall prediction skill in winter (cf. Fig. 14a and Fig. S5a) and loses a significant amount of skill in predicting internal variability in spring at 1- to 3-month lead time (cf. Fig. 14b and Fig. S6b). It confirms that retaining more modes in the model helps capture more internal variability and trends’ contribution in SIE predictions.
5. Summary and discussion
In this study, we developed a reduced-dimension stochastic model—namely, a linear Markov model—to predict Arctic SIC at the seasonal time scale. The model is constructed in the multivariate (MEOF) space with variables from the atmosphere, ocean, and ice field so that it captures the covariability of the Arctic climate system. Because of the complexity of Arctic geography and regional–seasonal characteristics of atmosphere–ocean–sea ice interaction, the determinants of the model’s variables and its number of MEOF modes are not unique. However, the differences in prediction skill from using different numbers of variables and MEOF modes are relatively small compared to the Markov model for Antarctic sea ice (Chen and Yuan 2004). Based on a trial-and-error method, we chose 3 variables (SIC, SST, and SAT) with 11 MEOF modes as the final model configuration, which was used in this study to assess model performance and make forward forecast. The results from 3 variables and 11 modes are similar to the results derived from 6 variables and 11 modes. The final choice of model configuration better fits the summer and fall prediction but may be overfitting for winter and spring prediction.
The linear Markov model shows good skill in predicting signs of Arctic sea ice concentration anomaly a few seasons in advance. The predictability strongly depends on seasons and locations and is especially high within the Arctic basin during summer and fall. The model AC stays above 0.6 even at a 9-month lead for these seasons in the Chukchi, Beaufort, Eastern Siberian, Kara, and Barents Seas. The winter skill is 0.5 for up to 9-month-lead forecasts in the marginal and peripheral seas of the Arctic, such as the Barents Sea, Greenland Sea, and Baffin Bay. The pattern of prediction skill in spring is very similar to that in winter but with lower correlations. In general, the forecast skill is better in the Atlantic sector of the Arctic than in the Pacific sector and better in summer and fall than in winter and spring. The lowest prediction skill occurs in the Bering Sea. On the pan-Arctic average, the model’s skill is superior to the skill derived from anomaly persistence, damped anomaly persistence, and climatology predictions. The long-term trends in the Markov model predictions contribute largely to the SIC summer and fall forecast skill in the Chukchi, East Siberian, and Barents Seas as well as winter forecast skill in the peripheral seas of the Atlantic sector. However, the model retains skillful 3-month-lead summer and fall predictions in the central Arctic and 6-month-lead fall predictions in the Canadian basin after the trends were removed. Moreover, the model’s skill is consistently higher than the persistence skill after all trends were removed from summer and fall predictions, revealing the model’s ability to capture some predictable SIC internal variability. The success of the model is attributed to the dominance of several distinct modes in the coupled atmosphere–ocean–sea ice system and to the model’s ability to pick up these modes. However, the relationships between sea ice and the ocean as well as atmosphere likely vary in different climate regimes (Holland and Stroeve 2011). A statistical model would not be at its best performance for a long period, particularly for the rapidly declining arctic sea ice. Periodically retraining the Markov model is needed to keep its skill updated.
In addition, the SIE prediction can be derived from predicted SIC. Our results show that SIE is most predictable in September and in winter and least predictable in spring (Fig. 14). Contrasting to the trend’s contribution to SIC prediction skill, the calculated SIE predictions cannot capture the rapid decline trends in observed SIE but are able to capture the internal variability. Removing monthly trends from SIE prediction improves its prediction skill. The SIE predictions without linear trends are superior to the detrended persistence for most months and lead times. Moreover, the SIE prediction by the Markov model consistently outperforms the climatology prediction in terms of RMSE except in spring, regardless of whether one considers trends. It demonstrates that the model predicts Arctic SIE’s internal variability better than its long-term trends. Predicting September SIE is given as an example to illustrate the model’s forecast skill for SIE. Although the model has a relatively high skill in predicting the Arctic SIE, it systematically overestimates the September SIE. This model bias can be corrected by two methods: a constant bias correction by removing the mean bias over the last 34 yr and a linear regression bias correction by fitting model errors as a linear function of initial SIE anomalies and then removing the regressed errors. The regression bias correction can significantly improve the skill of 2-month-lead forecasts and provide better prediction than the constant bias correction (Table 1 and Fig. 12). Because the Markov model for SIC is built in the multivariate space, it might underestimate the SIE trend since sea ice likely has the fastest declining trend in the system. The regression bias correction can rectify the trend quite well at 2-month lead time. On the other hand, the 4-month-lead forecast cannot benefit from the regression bias correction since prediction errors have a near-zero slope in the regression (Fig. 12).
This study shows that the SIC within the Arctic basin can be predicted seasons in advance by a reduced-dimension statistical model with reasonable skill. The Markov model effectively utilizes the linear relationships between SIC and oceanic and atmospheric variability. Future improvement of the Arctic sea ice prediction by statistical models likely relies on the exploration of the nonlinear component existing in the Arctic climate system, or non-Markovian processes, as well as on more sophisticated bias corrections than we presented in this study. Moreover, the low prediction skill in the Pacific sector of the Arctic may not be improved in a pan-Arctic model. One of reasons could be that the dominant climate variability in the northern mid–high latitudes occurs mostly in the Atlantic sector of the Arctic and subarctic, which dictates leading MEOF modes. The unique characteristics of atmosphere–ocean–sea ice coupled relationships in the Pacific sector may not be included in the leading MEOF decompositions of the pan-Arctic climate system and thus are not correctly represented in the model. On the other hand, the Pacific sector of the Arctic may need a different set of variables to maximize the model’s predictability. A separate study shows that the sea ice prediction in the Bering Sea and the Sea of Okhotsk can be significantly improved in a regional linear Markov model that focuses on the important coupled atmosphere–ocean–sea ice relationships in these regions (Li et al. 2016).
We gratefully acknowledge support from the Office of Naval Research Grant N00014-12-1-0911. Anonymous reviewers’ comments helped us greatly improve the quality of the study.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-15-0858.s1.
Lamont-Doherty Earth Observatory Contribution Number 8062.