Because sea ice thickness is known to influence future patterns of sea ice concentration, it is likely that an improved initialization of sea ice thickness in a coupled ocean–atmosphere model would improve Arctic sea ice cover forecasts. Here, two sea ice thickness datasets as possible candidates for forecast initialization were investigated: the Climate Forecast System Reanalysis (CFSR) and the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS). Using Ice, Cloud, and Land Elevation Satellite (ICESat) data, it was shown that the PIOMAS dataset had a more realistic representation of sea ice thickness than CFSR. Subsequently, both March CFSR and PIOMAS sea ice thicknesses were used to initialize hindcasts using the Climate Forecast System, version 2 (CFSv2), model. A second set of model runs was also done in which the original model physics were modified to more physically reasonable settings—namely, increasing the number of marine stratus clouds in the Arctic region and enabling realistic representation of the ice–ocean heat flux. Hindcasts were evaluated using sea ice concentration observations from the National Aeronautics and Space Administration (NASA) Team and Bootstrap algorithms. Results show that using PIOMAS initial sea ice thickness in addition to the physics modifications yielded significant improvement in the prediction of September Arctic sea ice extent along with increased interannual predictive skill. Significant local improvements in sea ice concentration were also seen in distinct regions for the change to PIOMAS initial thickness or the physics adjustments, with the most improvement occurring when these changes were applied concurrently.
According to the Fifth Assessment Report from the Intergovernmental Panel on Climate Change (IPCC), annual Arctic sea ice extent (SIE) is very likely (90%–100% confident) to have decreased at a rate of 0.45 to 0.51 million km2 decade−1 during the 1979–2012 period (Vaughan et al. 2013), leading to projections of a summer ice free Arctic by the 2030s (Wang and Overland 2012). Sea ice loss can be attributed to both anthropogenic influences and natural variability (Kay et al. 2011; Swart et al. 2015). On a seasonal time scale, accurate Arctic sea ice prediction is important for oil and shipping interests, wildlife protection, and ecosystems management.
Studies support that the Arctic sea ice cover is potentially predictable up to several years in advance (Blanchard-Wrigglesworth et al. 2011; Day et al. 2014a; Tietsche et al. 2014). However, assessments only show actual predictive skill for sea ice out to only a few months (Chevallier et al. 2013; Merryfield et al. 2013; Sigmond et al. 2013; Wang et al. 2013). The difference between potential predictability and actual skill raises interesting questions for possible causes—for example, errors in the initialization of sea ice thickness, which has been labeled as a good potential predictor of the sea ice cover (Lindsay et al. 2008). Wang et al. (2013) showed that changes in sea ice thickness during the spring months had a large impact on forecasts of September SIE using the Climate Forecast System, version 2 (CFSv2), with thickness increases resulting in higher SIE forecasts and vice versa, which agreed with observations. Low March 2007 sea ice thickness anomalies were cited by Kauker et al. (2009) as a factor leading to the record low sea ice cover in the following September. Day et al. (2014a,b) showed that sea ice thickness information is crucial for prediction of sea ice cover up to eight months in advance, and Chevallier and Salas-Melia (2012) demonstrated the area covered by thick ice above a certain threshold in March can be used to predict September sea ice area (SIA). Msadek et al. (2014) also showed skillful predictions of September SIE can be obtained from initializations as early as March using improved atmospheric initial conditions and highlighted the potential benefits of initializing improved sea ice thickness data.
Therefore, given that sea ice thickness is a potentially important predictor to forecast September sea ice cover, it is important to investigate the uncertainties in the initialization of sea ice thickness in the current forecast systems and to what extent improvements in the initial condition sea ice thickness dataset used in models could yield an improvement in September sea ice forecasts, thereby narrowing the gap between prediction skill and potential predictability. In this study, we analyze the sea ice thickness from two data assimilation systems and demonstrate the dependence of sea ice forecast skill on the accuracy of initial sea ice thickness based on hindcast experiments.
a. Sea ice thickness datasets
Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) sea ice thickness data (Zhang and Rothrock 2003) from the Polar Science Center at the University of Washington were used, in addition to Climate Forecast System Reanalysis (CFSR; Saha et al. 2010) sea ice thickness data from the National Centers for Environmental Prediction (NCEP). Sea ice components in both PIOMAS and CFSR assimilate satellite measurements of sea ice concentration. Sea ice thickness is derived from internal dynamics and thermodynamics without assimilating observed thickness information. CFSR uses the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (MOM4) coupled with the GFDL Sea Ice Simulator. PIOMAS is run with the Parallel Ocean Program coupled with a thickness and enthalpy distribution sea ice model. CFSR ocean and ice models are coupled with an atmospheric model and data assimilation system while PIOMAS ocean and ice models are forced with the NCEP–National Center for Atmospheric Research (NCAR) reanalysis. CFSR provides initial conditions for the NCEP operational seasonal climate prediction system, CFSv2 (Wang et al. 2013; Saha et al. 2014), and PIOMAS has been used to initialize sea ice forecasts for other models (Zhang et al. 2008; Blanchard-Wrigglesworth et al. 2015, manuscript submitted to Geophys. Res. Lett.). An assessment of the difference between CFSR and PIOMAS and its impact on sea ice prediction will help understand the source of errors in the prediction.
Figure 1 compares sea ice thickness and integrated sea ice volume from CFSR, PIOMAS, and the Ice, Cloud, and Land Elevation Satellite (ICESat; Schutz et al. 2005). ICESat data used are described in Kwok et al. (2009) and available for download at http://rkwok.jpl.nasa.gov/icesat. Better agreement is seen between the spatial map of ICESat sea ice thickness (Fig. 1c) and the spatial map of PIOMAS thickness (Fig. 1b) than the spatial map of CFSR thickness (Fig. 1a). Schweiger et al. (2011) suggested that ICESat thicknesses may be too high off the northern coast of Greenland (Fig. 1c) as PIOMAS data were found to be more comparable with in situ observations in this region. Laxon et al. (2013) also established that the representation of sea ice volume from PIOMAS is realistic based on comparisons with ICESat and CryoSat-2. This is also addressed in Fig. 1d, which shows the mean February and March sea ice volume from PIOMAS and CFSR for the years 1982–2013 with the ICESat data plotted at the appropriate times. It is apparent that the ICESat mean state is more in line with PIOMAS than with CFSR. With this information, the use of PIOMAS sea ice thickness data as initial conditions with goals of improving CFSv2 forecasts of sea ice cover is justified.
b. Model setup
To assess the influence of sea ice thickness on subsequent forecasts, a set of CFSv2 hindcasts were initialized in March 2009–13 and run for nine full target months (April–December) with both CFSR and PIOMAS initial sea ice thickness thus covering the sea ice minimum each year. For each year, five ensembles were run (initialized 0000 UTC 8–12 March). The length of nine target months follows the setting for NCEP operational seasonal prediction with CFSv2. The choice of March as initial months allowed an analysis for both melting in summer and refreezing after fall. First, we focused on the assessment of the evolution of the mean seasonal cycle across the five hindcast years (2009–13), in particular the representation of the September sea ice minimum. While a period of five years allowed for an evaluation in the prediction of seasonal cycle as well as an assessment of techniques for initializing forecasts with PIOMAS sea ice thickness and for generating an ensemble of forecasts, it was too short for a comprehensive assessment of the model performance in predicting interannual variability. Therefore, to analyze year-to-year prediction, select hindcasts were extended to cover 10 years (2005–14, section 3c).
CFSv2 is a coupled ocean–atmosphere model, with the components being the NCEP Global Forecast System (GFS) and GFDL MOM4 for the atmosphere and ocean, respectively (Saha et al. 2014). Control runs were set up with CFSR initial sea ice thickness (CFSv2CFSR). PIOMAS data were interpolated to CFSR’s Arctic grid spacing (as in CFSv2). In addition, PIOMAS ice thickness categories were converted from their native 12-category configuration (Zhang and Rothrock 2001) to the 5-category configuration used in CFSR (Saha et al. 2010). Following these interpolations, CFSv2 runs were initialized with the PIOMAS sea ice thickness (CFSv2PIOM) but with all other fields from CFSR.
A second set of hindcasts was run with a modification to the model physics (CFSv2CFSRp and CFSv2PIOMp). One major issue in the previous CFS, version 1 (CFSv1), was the lack of marine stratus clouds, which caused large sea surface temperature (SST) warm biases off the tropical west coasts of continents and mid- to high latitudes of the summer hemisphere (Moorthi et al. 2010). For the CFSR and CFSv2, a parameterization scheme was developed to improve the simulation of maritime clouds. The modified scheme includes two changes to the original atmospheric model physics: 1) limit the shallow convection top to be below the low-level inversion when the condition for cloud-top entrainment instability is not satisfied and 2) set the background vertical diffusion to zero above low-level inversions. The combination of these two modifications leads to an improved prediction of marine stratus (Moorthi et al. 2010). This scheme has been used in the CFSR (Saha et al. 2010).
During the testing and development phase of CFSv2, the forecast system counterpart of CFSR, this scheme was disabled because its use together with the easterly wind bias in the tropical Pacific in the atmospheric component caused a large SST cold bias and greatly weakened the predicted El Niño–Southern Oscillation (ENSO) variability (Saha et al. 2014). This step, resulting in essentially turning off the generation of stratus clouds, was justified as the performance in ENSO prediction was the dominant concern, while less attention was paid to the skill of sea ice prediction. Disabling of the stratus cloud scheme in the CFSv2 was for a pragmatic reason: to maintain reasonable tropical SST variability rather than being physically based. However, this also resulted in a large SST warm bias in the high latitudes during warm seasons due to excessive surface downward solar radiation, causing too-rapid sea ice melting. To compensate for this, an artificial upper limit for the bottom heat flux from ocean water into sea ice was imposed—another unphysical treatment in the CFSv2 model. In this study, physics modifications refer to 1) “reenabling” the scheme of Moorthi et al. (2010) to improve marine stratus clouds and 2) “removing” the constraint for water–ice heat flux to revert the model configuration to what it should physically be. Comparison between these two hindcast sets allows an analysis of the dependence of the impact of sea ice thickness initialization on these model physics settings. A list of model runs integrated in this study is provided in Table 1.
Model hindcasts were evaluated against observations from the National Aeronautics and Space Administration (NASA), specifically the NASA Team (Cavalieri et al. 2014) and NASA Bootstrap (Comiso 2014) sea ice concentrations, which were interpolated to match the gridcell spacing of the CFSv2 output. Both datasets are dependent on the same set of satellite observations, but differences arise in the interpretation of meltwater processes. As explained in Notz (2014), summer sea ice has higher values of sea ice concentration in the NASA Bootstrap dataset than in the NASA Team dataset. Although both algorithms treat sea ice covered by surface meltwater as open water, the Bootstrap algorithm compensates for this bias more than the NASA Team algorithm. The NASA observations have a region close to the pole that cannot be observed owing to the orbit inclination of the satellites. This region, known as the “polar hole” with an area of approximately 0.28 × 106 km2, is removed from all datasets for a concise evaluation.
The assessment of the model output was done using hemispheric scale total sea ice cover in addition to gridcell sea ice concentration using changes in absolute error of the different model configurations relative to the base scenario (absolute error of CFSv2CFSR with respect to NASA observations). Significant improvements in error at 99% confidence were determined using a one-tailed t test. Total sea ice cover is commonly represented by SIE and SIA. Following the IPCC report (Vaughan et al. 2013), SIE is defined as the region encompassed by the edge of sea ice, represented by a concentration of at least 15%. Because of this, regions with small leads in the sea ice with a concentration greater than the 15% threshold are still counted as part of the SIE. SIA is different, as these leads are not included. For this reason SIE is always greater than SIA. In this analysis, Arctic SIE was computed by taking the cumulative sum across the Northern Hemisphere of the area of each grid cell with a sea ice concentration at or above 15%. SIA was calculated by multiplying the sea ice concentration in each grid cell by the area of the respective grid cell and taking a Northern Hemisphere sum of the results for grid cells with a sea ice concentration at or greater than 15%. For each hindcast year, SIE and SIA were calculated individually for each ensemble member and then averaged to yield a mean SIE and SIA value for that year. Finally, rather than taking the hemispheric total, sea ice concentration values at individual grid cells were also examined to better assess regional patterns.
a. Integrated SIE and SIA
We first evaluated the seasonal cycle of SIE and SIA from the experiments. Figure 2 presents the monthly mean forecast SIE and SIA (excluding the polar hole region). Figure 3 shows the month to month change of each variable to illustrate how well the seasonal cycle is captured. Absolute differences in SIE and SIA relative to NASA observations are shown in Fig. 4.
Focusing on total SIE, it is evident that either modifying the model physics or changing the initial thickness to PIOMAS resulted in improved forecasts of the sea ice minimum relative to CFSv2CFSR. Prediction of mean September 2009–13 sea ice minimum in the base CFSv2CFSR run yielded a SIE value of 7.96 × 106 km2. Actual observations from NASA Team and NASA Bootstrap were 5.04 × 106 and 4.97 × 106 km2, respectively. When only modifying the physics and still using CFSR thickness initialization, the predicted mean September SIE dropped to 6.35 × 106 km2. Switching to PIOMAS initialized thickness the predicted SIE was 6.31 × 106 km2 without the physics change and 4.72 × 106 km2 with the physics change. Clearly, based on Fig. 4, using PIOMAS in conjunction with the physics change produced the best prediction of September mean SIE with significant improvements at 99% confidence relative to CFSv2CFSR for both comparisons with NASA Team and NASA Bootstrap. For the NASA Team comparison the improvement was significant throughout the July–October period (Fig. 4a) and September for NASA Bootstrap (Fig. 4b). The increase in absolute error in August SIE in the CFSv2PIOMp runs can be attributed to the low prediction of SIE relative to the observations (Fig. 2c).
SIA was more difficult to evaluate given the differences that exist within the observational datasets, with NASA Bootstrap consistently producing a higher SIA than NASA Team because it has larger concentration values (Notz 2014). There are also seasonal differences in the minimum in SIA with the runs with physics changes producing the minimum in August and the runs with no physics changes and the observations showing the minimum in September (Figs. 2d,f). For this reason, we focus on the average August–September value. August–September average SIA output from CFSv2PIOM and CFSv2PIOMp were 4.40 × 106 and 3.78 × 106 km2, respectively. This compares better to the NASA observations (3.48 × 106 km2 for Team and 4.34 × 106 km2 for Bootstrap) than the runs with CFSR initialized thickness (5.68 × 106 km2 for CFSv2CFSR and 5.11 × 106 km2 for CFSv2CFSRp). Absolute error changes with respect to NASA Team are significantly decreased for CFSv2PIOMp for July–September but only in September for CFSv2PIOM (Fig. 4c). Relative to NASA Bootstrap, it is apparent from Fig. 4d that CFSv2PIOMp shows a continual improvement while the other configurations have an increase in error between July and November. The improvement in CFSv2PIOMp is particularly significant in November.
As mentioned before, part of the issue also arises from how the different model configurations handle the seasonal cycle of sea ice. As shown in Figs. 3c and 3f, the acceleration of ice formation (represented by the slope of the plotted lines) in the August–October period better agrees with the NASA observations in runs which used both changes in physics and PIOMAS initial thickness than the other runs. Upon closer examination, it is evident that the PIOMAS initial sea ice thickness conditions allow for a better representation of the SIE and SIA tendency in August–September (melt; Figs. 3b,e) and the physics changes show higher skill in September–October (freeze; Figs. 3a,d). Without the PIOMAS initial sea ice thickness the melting rate is too slow, and without the physics changes, the refreeze is too slow compared to observations. The optimal configuration is when both physics and initial thicknesses are changed (Figs. 3c,f). The only caveat is that the actual change in SIA is positive in September for the runs with physics changes and negative in the observations (Figs. 3d,f), which leads to the discrepancy in the modeled minimum as explained in the preceding paragraph. The next section will focus on sea ice concentration on the grid cell scale to address the SIA prediction uncertainties.
b. Regional sea ice concentration
For local characteristics among the experiments, we first looked at prediction around the Bering Strait, where the sea ice seasonal cycle is largely controlled by thermodynamic processes (Bitz et al. 2005) and thus strongly affected by initial sea ice thickness and surface heat fluxes. Figure 5 shows zonal mean sea ice concentration from April to December averaged for 2009–13 between 170° and 200°E. The NASA Team and NASA Bootstrap data were essentially the same with a strong seasonal cycle (Figs. 5e,f). Taking the location of 15%–30% concentration to measure the sea ice evolution, the observed sea ice retreated northward from near 65°N in June to 76°N in August and September and then expanded southward reaching 62°N in December. CFSv2CFSR produced weak seasonal cycle with 15%–30% concentration reaching slightly to the north of 70°N in September and returning to around 66°N in December (Fig. 5a). Improvements were seen with the modifications to the model physics or with the use of PIOMAS sea ice initial conditions. In September, the sea ice retreated to 73°N in CFSv2CFSRp (Fig. 5b) and 75°N in CFSv2PIOM (Fig. 5c) and CFSv2PIOMp (Fig. 5d), suggesting the use of PIOMAS initial thickness conditions created more improvement than just the change to model physics in this region since the NASA observations showed retreat up to 76°N. It is noted that the delayed autumn sea ice formation in CFSv2CFSR is not improved in the other experiments, suggesting that this bias in CFSv2 is not related to the initialization of sea ice thickness or the physics tested in this study.
Next, we compared spatial patterns of September sea ice concentration throughout the entire Arctic. Figure 6 shows the mean September sea ice concentration across the region for the four model configurations and the NASA observations, with the runs with physics changes and the NASA Bootstrap observations having the highest values supporting the SIA analysis in the previous section. As apparent in Fig. 7, there are three distinct regions that showed significant improvements in absolute error relative to NASA observations—namely, the North Atlantic, the central Arctic, and the northern coasts of Alaska and Russia. The regions with changes are a function of the model configuration used and the observational dataset compared to. For example, applying physics changes significantly improved central Arctic sea ice forecasts relative to NASA Bootstrap (Figs. 7g,h). Adjusting the physics improved the hindcasts in the North Atlantic relative to both observations (Figs. 7c,d,g,h). Changing the thickness initiation improved hindcasts off the northern Alaska and Russian coasts for both CFSv2PIOM and CFSv2PIOMp (Figs. 7b,f,d,h), which agrees with the results presented in Fig. 5. When both the physics were adjusted and PIOMAS thickness was used to initialize the model, the significant improvements were seen in both the North Atlantic and off the Alaskan–Russian coast, yielding this as the optimal configuration. Central Arctic improvements are only dependent on the observed dataset used to evaluate but not the initial thickness or physics settings. The total area of the SIE region with significant improvement in absolute error relative to CFSv2CFSR is presented in Fig. 7i and Fig. 7j using NASA Team and NASA Bootstrap observations, respectively, as a function of month, and it is apparent that using PIOMAS sea ice thickness with adjusted physics yields the largest region of improvements in hindcasts throughout most of the period.
c. Interannual variability
Finally, it is pertinent to look into interannual variability. For this analysis, the hindcasts for CFSv2CFSR and CFSv2PIOMp were extended to cover 2005–14. Only the hindcasts for these two configurations were extended for this part of the analysis as they represent the lowest and highest skill, respectively, in terms of the mean state prediction as shown in the previous sections. Here we used September SIE and SIA from individual ensemble members and looked into the mean and spread for individual years rather than a seasonal average. Figure 8 shows a time series for the two CFSv2 model configurations as well as the NASA observations. It is apparent that for both SIE and SIA that CFSv2PIOMp agrees better with the NASA observations than CFSv2CFSR for each individual year. Tables 2 and 3 provide statistical parameters for SIE and SIA, respectively.
The overall linear trend in SIE and SIA was more accurately captured in CFSv2PIOMp than in CFSv2CFSR relative to the NASA observations as illustrated by the dotted lines in Fig. 8 and the first columns of Tables 2 and 3. The trend is near neutral in SIE and SIA from CFSv2CFSR while CFSvPIOMp and both sets of NASA observations show trends in SIE and SIA between −0.06 × 106 and −0.09 × 106 km2 yr−1. Higher correlation coefficients and lower root-mean-square errors relative to the NASA observations were found when the CFSv2PIOMp configuration was used than from using CFSv2CFSR. Even when linear trends were removed from all datasets, this holds true, indicating some degree of improved prediction of sea ice interannual variability was achieved using CFSv2PIOMp. Despite the marked statistical improvement in the CFSv2PIOMp hindcasts, the ensemble spread SIE hindcasts only encompassed the NASA observations in 4 out of 10 years. This was not as clear for SIA because of the larger variability between NASA Team and Bootstrap observations.
4. Discussion and conclusions
Hindcasts were conducted with the CFSv2 model initialized in March using different sea ice thickness datasets and changed model physics parameters to assess improvements in sea ice prediction. While previous studies (Wang et al. 2013; Msadek et al. 2014) have pointed out the difficulty in using a more realistic sea ice thickness dataset due to the sparseness of observations, PIOMAS output, which covers the entire Arctic region, was found to be observationally consistent and therefore determined as a feasible dataset to use in this study. Despite the short hindcast period, significant results at 99% confidence were still able to be found. The results show that hindcasts of SIE were improved when PIOMAS sea ice thickness was used as the initial conditions for CFSv2. In addition, the hindcasts were further improved when modifications to the physics were made to restore the model to physically reasonable settings, and the improvements were significant at 99% confidence for the month of September. SIA proved harder to evaluate owing to differences in the observation datasets and NASA Team and NASA Bootstrap satellite sea ice concentration retrievals. However, the overall tendency was best modeled in the runs that both adjusted the physics and used PIOMAS initial thickness conditions. It is also worth noting that the unrealistic positive trend in early spring sea ice volume in CFSR (Fig. 1d) corresponded to a neutral trend in predicted September SIE and SIA (Fig. 8) while the downward sea ice volume trend seen in PIOMAS during February and March translated into a more observationally consistent downward trend in predicted SIE and SIA.
Seasonal evolution of sea ice concentration in local regions showed similar behavior among the hindcasts. In particular, the original CFSv2 initialized from CFSR produced a weak seasonal cycle of sea ice concentration around the Bering Sea and Chukchi Sea. The change in physics and the use of PIOMAS sea ice thickness as the initial conditions resulted in an improved sea ice seasonal cycle. A closer look into September sea ice concentration across the Arctic identified relative impacts from changes in initialization compared to those from changes in model physics. Specifically, physics changes significantly improved sea ice prediction within the North Atlantic, initial thickness changes improved the forecast off the northern Alaska and Russian coasts, and both physics and initial thickness changes resulted in improvements in both aforementioned regions highlighting this as the optimal configuration. The North Atlantic and coastal Alaska/Russia regions are along the edge of the sea ice, which shows that the model enhancements appear proficient at better predicting the boundary between sea ice and open ocean, but at this time because of the uncertainties discussed with SIA, it cannot be said with confidence that there are significant improvements in forecasting the exact amount of sea ice that is present within the ice covered region.
The results presented here call for accurate initialization of sea ice thickness as well as realistic physical treatments in the coupled system for more accurate seasonal prediction of the seasonal cycle of sea ice as they can benefit both the prediction of the mean state and interannual variability. Improved knowledge of the location of sea ice edge will greatly benefit many operations based in the Arctic. Future work will need to assess the dependence of the improvement on initial months. For example, it is possible that an accurate initial condition of sea ice thickness may benefit hindcasts and forecasts of September sea ice cover initialized from June more than that from March (Chevallier and Salas-Melia 2012; Day et al. 2014a). Work is also planned to address the predictability of the first calendar year/day of sea ice melt and freeze using this new initialization, particularly over the Alaska region where the initial sea ice thickness change appeared to have the largest positive results.
Support for JZ was provided by ONR Grant N00014-12-1-0112. PIOMAS data used in this study are available at ftp://pscftp.apl.washington.edu/zhang/IDAO/retrospection and NASA Team and Bootstrap sea ice concentrations can be downloaded at ftp://sidads.colorado.edu/pub/DATASETS. We thank Meng-Pai Hung (CPC) for processing ICESat data used in this study.