Here sea ice concentration derived from the Special Sensor Microwave Imager/Sounder and thickness derived from the Soil Moisture and Ocean Salinity and CryoSat-2 satellites are assimilated in the National Centers for Environmental Prediction Climate Forecast System using a localized error subspace transform ensemble Kalman filter (LESTKF). Three ensemble-based hindcasts are conducted to examine impacts of the assimilation on Arctic sea ice prediction, including CTL (without any assimilation), LESTKF-1 (with initial sea ice assimilation only), and LESTKF-E5 (with every 5-day sea ice assimilation). Assessment with the assimilated satellite products and independent sea ice thickness datasets shows that assimilating sea ice concentration and thickness leads to improved Arctic sea ice prediction. LESTKF-1 improves sea ice forecast initially. The initial improvement gradually diminishes after ~3-week integration for sea ice extent but remains quite steady through the integration for sea ice thickness. Large biases in both the ice extent and thickness in CTL are remarkably reduced through the hindcast in LESTKF-E5. Additional numerical experiments suggest that the hindcast with sea ice thickness assimilation dramatically reduces systematic bias in the predicted ice thickness compared with sea ice concentration assimilation only or without any assimilation, which also benefits the prediction of sea ice extent and concentration due to their covariability. Hence, the corrected state of sea ice thickness would aid in the forecast procedure. Increasing the number of ensemble members or extending the integration period to generate estimates of initial model states and uncertainties seems to have small impacts on sea ice prediction relative to LESTKF-E5.
Arctic sea ice extent and thickness have experienced dramatic change in the past few decades. Sea ice extent has declined for all months since the late 1970s (e.g., Comiso 2012; Cavalieri and Parkinson 2012); that is, September ice extent has declined 13.3% decade−1 during 1979–2016, which is underestimated by most of the global climate models that participated in phase 5 of the Coupled Model Intercomparison Project (CMIP5) (Stroeve et al. 2012). Accompanying the rapid decline of the ice extent has been a thinning of the ice pack; that is, compared to the submarine data during 1958–76, the Ice, Cloud, and Land Elevation Satellite (ICESat) data during 2003–08 show that the mean ice thickness has decreased by ~1.6 m (e.g., Kwok and Rothrock 2009). This is largely a result of thinner first-year ice replacing thicker multiyear ice (e.g., Maslanik et al. 2011). Rapid change of Arctic sea ice has increased interannual variability in the ice extent, particularly from summer to autumn (Serreze and Stroeve 2015). The decreasing Arctic sea ice has also coincided with a period of more frequent extreme weather events in midlatitudes of the Northern Hemisphere (i.e., cold and snowy winters over parts of Eurasia and North America) (e.g., Liu et al. 2012; Cohen et al. 2012, 2014; Francis and Vavrus 2012; Mori et al. 2014). Analyses of climate model projections suggested that the summer Arctic would reach an ice-free state in the middle of this century under high-emission scenarios (e.g., Wang and Overland 2012; Liu et al. 2013; Notz and Stroeve 2016), which offers faster routes for commercial shipments between Europe/North America and Asia (Smith and Stephenson 2013; Melia et al. 2016). Rapid change of Arctic sea ice has captured attention and posed significant challenges to a wide range of stakeholders, including maritime safety and security, resource management and development, coastal communities, climate change researchers, politicians, and a growing segment of the general public (e.g., Stroeve et al. 2014). Hence there is a rising demand for Arctic sea ice prediction (e.g., Eicken 2013).
Seasonal prediction of Arctic sea ice has been primarily produced with statistical methods in the past (see http://www.arcus.org/sipn). Some predictions use ensemble simulations from coupled sea ice–ocean models with prescribed atmospheric forcings (e.g., Kauker et al. 2009; Zhang et al. 2008). Recently, a few operational centers have implemented a sea ice model component in their climate forecast systems for seasonal climate prediction, for example, NCEP Climate Forecast System, version 2 (CFSv2; Saha et al. 2014). Sea ice prediction is particularly challenging in the context of fully coupled climate forecast system. In a fully coupled climate forecast model, changes in sea ice concentration and thickness regionally or globally would influence atmospheric and oceanic conditions, which in turn affect sea ice distribution.
As suggested by sea ice prediction network (Stroeve et al. 2014), large biases exist in the predicted sea ice extent and thickness. One of the primary sources of uncertainty in seasonal sea ice prediction is poorly known initial sea ice conditions. To reduce such uncertainty, reliable sea ice observations and a feasible data assimilation system are required to produce reasonable initial conditions.
Because of the extensive spatial and temporal coverage of sea ice concentration, several studies have explored assimilating satellite-derived ice concentration into stand-alone sea ice models, coupled sea ice–ocean models, and fully coupled climate models using different assimilation methods: ensemble Kalman filter by Lisæter et al. (2003), observational nudging scheme by Lindsay and Zhang (2006), optimal interpolation by Stark et al. (2008), three-dimensional variational method by Caya et al. (2010), and local ensemble singular evolutive interpolated Kalman filter (SEIK) by Yang et al. (2014, 2015). These studies demonstrated that the assimilation of observed sea ice concentration in the model improves the simulation of sea ice cover. Some of these studies also noticed that improvement in the simulation of sea ice thickness is very small.
To better predict sea ice, accurate sea ice initialization requires not only sea ice concentration but also variables (i.e., sea ice thickness) that influence surface energy fluxes and, thereby, ocean–atmosphere interactions. Using the June 2015 sea ice prediction from the NCEP CFSv2 as an example, the forecasted September Arctic sea ice extent based on an ensemble mean of the NCEP CFSv2 was 5.7 × 106 km2 (an estimated error of ±0.47 × 106 km2). This forecast removed the model’s systematic bias based on a statistical regression between what the NCEP CFSv2 predicts and what is observed. By contrast, the forecasted September ice extent from the ensemble mean of the NCEP CFSv2 with revised May 2015 initial conditions was 4.6 × 106 km2 (an estimated error of ±0.44 × 106 km2). The initial condition was modified from the real time of the NCEP CFSv2 each day by thinning the ice, since it is known that the NCEP CFSv2 is biased toward sea ice being too thick. The latter method is based on thinning the ice and seeing whether the extent is thicker than a critical limit. The large discrepancy between two sea ice seasonal predictions suggests a strong impact of sea ice thickness assimilation (http://www.arcus.org/sipn/sea-ice-outlook/2015/june; Day et al. 2014a). At seasonal time scales, the preconditioning of spring sea ice thickness has been shown to be important for summer sea ice prediction (e.g., Wang et al. 2013; Day et al. 2014a,b). Higher predictability of sea ice thickness with respect to that of sea ice cover has also been found for at least 2 yr (e.g., Holland et al. 2011; Guemas et al. 2014).
Nevertheless, observing sea ice thickness is challenging owing to its spatial and temporal inconsistency (e.g., Kwok and Sulsky 2010). As a consequence, there are few assimilation studies using satellite-derived sea ice thickness, especially in fully coupled climate forecast models, although it has potential for reducing the bias and improving sea ice prediction as mentioned above. The purpose of this study is to assimilate recently computed and validated satellite-based Arctic sea ice thickness (as well as sea ice concentration) using a recently developed localized error subspace transform ensemble Kalman filter (LESTKF), which provides sea ice initial conditions for a fully coupled climate forecast system. Numerical experiments are then conducted to examine the impacts of assimilating sea ice thickness and concentration in combination as well as assimilating sea ice concentration only on the Arctic sea ice prediction in a fully coupled climate forecast system.
2. Model and data
The CFSv2 is a fully coupled atmosphere–land–ice–ocean model, which was made operational at the NCEP in 2011 (Saha et al. 2014). The atmospheric component is the NCEP Global Forecast System with a horizontal resolution of T126 (~1°) and 64 sigma–pressure hybrid layers vertically. The ocean component is the Geophysical Fluid Dynamic Laboratory (GFDL) Modular Ocean Model, version 4 (MOM4), with a longitudinal resolution of 0.5°, a latitudinal resolution of 0.25° between 10°S and 10°N, gradually increasing through the tropics until becoming fixed at 0.5° poleward of 30°S and 30°N, and 40 levels vertically (Griffies et al. 2004). A tripolar grid is utilized in MOM4 (Murray 1996), in which two poles are located at the landmasses of North America and northern Siberia, in addition to one at the South Pole. MOM4 includes a dynamic and thermodynamic sea ice model, the Sea Ice Simulator (SIS). The thermodynamics of the SIS is similar to Semtner’s thermodynamics with three layers vertically (one layer in snow and two layers in sea ice; Winton 2000). Five categories of sea ice thickness (0–0.1, 0.1–0.3, 0.3–0.7, 0.7–1.1, and >1.1 m) are used for each grid. A simple scheme transports ice between categories when the bounds of thickness in categories are exceeded because of thermodynamic and/or dynamic changes. The dynamics of the SIS uses the elastic–viscous–plastic rheology to calculate ice internal stress (Hunke and Dukowicz 1997).
The daily Arctic sea ice concentration is obtained from the National Snow and Ice Data Center, which is derived from the Special Sensor Microwave Imager/Sounder (SSMIS), processed by the bootstrap algorithm (Comiso 2000). Sea ice concentration used here is on the SSMIS polar stereographic grid with a spatial resolution of 25 km. Tonboe and Nielsen (2010) suggested that the average uncertainty of sea ice concentration is ~(10%–15%). Thus, a value of 0.15 is used to represent observational error, which is also commonly used as the threshold to define sea ice extent.
The observation of Arctic sea ice thickness is limited in space and time. Two satellite-derived sea ice thickness products are used here. One product is daily sea ice thickness derived from the Soil Moisture and Ocean Salinity (SMOS) brightness temperature using a single-layer emissivity model (Kaleschke et al. 2012; Tian-Kunze et al. 2014). The SMOS sea ice thickness, version 2, is used here (https://icdc.cen.uni-hamburg.de/thredds/catalog/ftpthredds/smos_sea_ice_thickness/v2/catalog.html), which is on the NSIDC polarstereographic grid with a spatial resolution of 12.5 km (Kaleschke et al. 2016). The uncertainties provided with SMOS ice thickness dataset are mainly from three sources: uncertainties of SMOS brightness temperature measurements, uncertainties in the auxiliary dataset, and assumptions made in the retrieval algorithm. The SMOS ice thickness has small uncertainty for thin ice but large underestimation for thick ice because of the saturation of SMOS brightness temperature with thickness (Kaleschke et al. 2012; Tian-Kunze et al. 2014). The averaged uncertainty of the SMOS ice thickness (<1 m) is ~0.7 m.
The other product is monthly sea ice thickness derived from freeboard measurements from the ESA’s CryoSat-2 satellite mission (e.g., Laxon et al. 2013). The CryoSat-2 sea ice thickness (obtained from http://data.seaiceportal.de) is on the NSIDC polar-stereographic grid with a spatial resolution of 25 km. The uncertainties of the CryoSat-2 ice thickness mainly come from freeboard and inadequate knowledge of snow depth and the radar interaction with snow (Ricker et al. 2014; Armitage and Ridout 2015). The CryoSat-2 random ice thickness uncertainty is up to 0.3 m considering the freeboard-to-thickness conversion, and its systematic ice thickness uncertainty is ~0.6 m (~1.2 m) for first-year (multiyear) ice (Ricker et al. 2014).
In this study, we use a simple combination of SMOS and CryoSat-2 sea ice thickness. Sea ice with thickness less than 1 m in SMOS is retained or otherwise replaced by CryoSat-2. As discussed above, the averaged uncertainty of the SMOS ice thickness (<1 m) is ~0.7 m, and the largest uncertainty for the CryoSat-2 ice thickness is ~1.5 m (~0.3 m for random and +~1.2 m for systematic). Here the largest uncertainty (1.5 m) is used for both the SMOS and CryoSat-2 sea ice thickness.
To assess our assimilation results, independent in situ sea ice thickness data are needed. Sea ice thickness measured from the NASA IceBridge mission is used for the evaluation. IceBridge is an airborne mission with the primary goal of bridging the gap between the end of ICESat in 2009 and the launch of ICESat-2 (Kurtz et al. 2013). Equipped with laser and radar instruments, IceBridge provides accurate measurements of sea ice thickness, particularly over level ice. The uncertainty of the measured ice thickness is ~0.4 m (Farrell et al. 2012), which is much smaller than that of satellite retrievals. The ice thickness from the NASA IceBridge sea ice freeboard, snow depth, and thickness quick look dataset is used here (https://nsidc.org/data/docs/daac/icebridge/evaluation_products/sea-ice-freeboard-snowdepth-thickness-quicklook-index.html), covering 14 March to 2 April in 2012, 21 March to 25 April in 2013, and 12 March to 28 April in 2014. IceBridge covers the Beaufort Sea, the Canadian Archipelago, and north of Greenland. Sea ice draft from sea ice mass balance buoy (IMB) is also used to assess the simulated ice thickness. Sea ice thickness is estimated from IMB by measuring the position of the ice bottom and the snow/ice surface using two sonar systems. The accuracy of the estimated ice thickness is ~0.1 m (Melling et al. 1995; Richter-Menge et al. 2006). The ice thickness from the ice mass balance buoy program is used here (http://imb-crrel-dartmouth.org/imb.crrel/buoysum.htm; Perovich et al. 2017). During the studying period, only one buoy (2013J) provides data outside the coverage of IceBridge, which initially deployed in the Laptev Sea and drifted toward the North Pole from 1 March to 15 April in 2013 (see Fig. 7 for its trajectory).
3. Assimilation method and numerical experiment
a. Assimilation method
An ensemble-based LESTKF is employed to assimilate satellite sea ice data and is implemented in the parallel data assimilation framework (PDAF; Nerger and Hiller 2013). Nerger et al. (2012a) showed that the error subspace transform Kalman filter (ESTKF) not only is a variant of SEIK (Pham 2001) but also has identical ensemble transformation as that of the ensemble transform Kalman filter (ETKF; Bishop et al. 2001). ESTKF provides consistent projections between the ensemble space and the error subspace as well as minimal ensemble transformation of the ensemble members, which outperforms SEIK and has lower computation coast (Nerger et al. 2012a). To reduce the sampling error of the estimated error covariance matrix, particularly long-range error covariance due to small ensembles, a regulated localization scheme is applied. This scheme conducts observation localization by weighting of the observation error covariance matrix through a localization function of variable width. For small errors, it avoids the widening of the effective localization length that may worsen the assimilation. Moreover, its additional computational cost is negligible relative to that of the analysis step (see details in Nerger et al. 2012b).
LESTKF consists of four steps: initialization, forecast, analysis, and resampling. Using the year 2012 as an example (one of the experiments in this study), for the initialization step, we generate initial ensembles by integrating the CFSv2 for one month (from 15 February to 15 March), which starts with an optimal analyzed initial condition generated by the NCEP Climate Forecast System (the data are available at https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/climate-forecast-system-version2-cfsv2#CFSv2 Operational Forecasts). For each day, the model state vectors of sea ice are stored sequentially to form the state matrix. The second-order exact sampling (Pham 2001) is then applied to the leading modes of empirical orthogonal functions (EOFs) of sea ice concentration in the Arctic to generate nine perturbations overlaid on the mean state of model trajectories. The nine ensembles provide an estimate of initial model states and uncertainties prior to the evolution of forecasts. Once the initial ensembles are generated, the forecast step begins, in which all ensembles are dynamically evolved with the CFSv2. At the end of the forecast step, the LESTKF filter analysis is conducted. For the analysis step, the aforementioned satellite-derived Arctic sea ice concentration and thickness are assimilated into the CFSv2, considering the observational errors and background uncertainty represented by the spread of model realizations. The analysis step updates the model mean state and covariance. After that, the ensemble resampling step is conducted, and all ensembles are updated to new states using the second-order exact sampling (Pham 2001) while preserving the ensemble mean and covariance.
b. Numerical experiments
Numerical experiments are performed using the CFSv2 to examine impacts of the assimilation of satellite-derived sea ice concentration and thickness on sea ice forecasts. Three sets of experiments are conducted. The first experiment is the control simulation without any assimilation (CTL), which is initialized with the state ensembles generated in the LESTKF initialization step. The second experiment is the same as CTL, except that the assimilation of sea ice concentration and thickness is performed only once at the beginning, and then the simulation is continued without further assimilation (LESTKF-1). The third experiment is the simulation by repeating the assimilation of sea ice concentration and thickness every 5 days (LESTKF-E5). In LESTKF-E5, the forecast, analysis, and resampling steps alternate until the end of the simulation. Three hindcasts are carried out from 7 March to 15 April for the years 2012, 2013, and 2014, respectively, given the availability of IceBridge observations.
As mentioned previously, the sea ice model component (SIS) in the CFSv2 uses five categories of sea ice thickness to describe sea ice thickness distribution. Thus, each category has a value of fractional ice area. The aggregated sea ice concentration of a grid cell is the sum of fractional ice areas for each category , (where ai is sea ice concentration in each thickness category and n = 1, …, 5). Similarly, the aggregated sea ice thickness of a grid cell is (where ai and hi are sea ice concentration and thickness for each thickness category and n = 1, …, 5). Satellite observations only provide the aggregated value of the ice concentration and thickness. After they are assimilated, we use a simple approach to assign the new aggregated ice concentration and thickness to the five categories as the initial conditions. That is, before the assimilation, the percentage of the ice concentration and thickness in each thickness category with respect to the aggregated value is calculated and stored. After the assimilation, the new aggregated ice concentration and thickness are reassigned to the five thickness categories based on the precalculated percentage. If the newly assigned ice thickness in one category is greater (less) than the corresponding ice thickness bound, then this category is shifted to the next (previous) category to meet the thickness limitation.
In the analysis step, a small influential cutoff distance of two model grid points (~100 km) is used in the regulated localization scheme. In addition, the statistical update of the state estimate under the assumption of Gaussian errors in the analysis step does not have a constraint on sea ice concentration and thickness. As a result, small negative and greater than 1 ice concentration can occur (i.e., at locations with small forecast concentration and no ice concentration observation). In such cases, small negative values are reset to 0, and the ice concentrations greater than 1 are reset to 1. For any grid point, if the ice thickness is positive but the ice concentration is 0, then the ice thickness is set to 0. If the ice concentration is positive but the ice thickness is 0, a scheme similar to Tietsche et al. (2012) is used, in which the ice thickness of 2 m is multiplied by the ice concentration. The sea surface temperature is not modified directly by the assimilation but is updated implicitly through sea ice and ocean coupling in the CFSv2. A forgetting factor of 0.98 is used to inflate the model error covariance after each data assimilation [see Nerger et al. (2012a) for details].
a. Sea ice extent and concentration
Figure 1 compares temporal evolution of the hindcast of the ensemble mean of the total Arctic sea ice extent (left panel) and the averaged root-mean-square error (RMSE) of sea ice concentration in each grid cell (right panel). Here the ice extent is calculated as the sum of the area of each grid cell having ice concentration greater than 15%. For consistency, we only calculate RMSE of sea ice concentration for the grid cell having ice concentration greater than 15%. At the initial stage, the simulated sea ice extent in CTL differs from the observation to varying degrees, very close to the SSMIS in 2012, much (relatively) less than the SSMIS in 2013 (2014). The deviation in the ice extent between the CTL and SSMIS becomes larger after several days of integration. Initially, the simulation with one-time sea ice assimilation (LESTKF-1) provides adjusted ice extent relatively closer to the observation compared to that of CTL (orange line vs red line). The simulated ice extent in LESTKF-1 is also relatively closer to the SSMIS than that of CTL through the model integration, although there is a tendency toward that of CTL after ~3-week integration. By contrast, LESTKF-E5 (every 5-day sea ice assimilation) obviously improves the ice extent simulation, in which the deviation from the SSMIS is greatly reduced through the hindcast (green line).
The RMSE of sea ice concentration (right panel in Fig. 1) shows that at the initial state as well as the first 2-week integration, large RMSE in CTL is obviously reduced as a result of the initial sea ice assimilation in LESTKF-1. After that, LESTKF-1 shows a similar level of RMSE as that of CTL (orange line vs red line). By contrast, LESTKF-E5 produces a zigzag RMSE time series (green line), having much (relatively) smaller RMSE relative to that of CTL once satellite sea ice data are assimilated (after ~3-day integration without further assimilation). Spatially, the above improved hindcast of sea ice extent and reduced RMSE are primarily achieved in the marginal ice zone in the North Atlantic and Pacific sectors (not shown). Here we examine the two-tailed statistical t test to determine whether differences among the above three hindcasts are significant. We consider that the time series of RMSE for the three hindcasts are independent random samples in normal distribution (Yang et al. 2015). Hereafter P1, P2, and P3 refer to a statistical p value between LESTKF-1 and CTL, between LESTKF-E5 and CTL, and between LESTKF-E5 and LESTKF-1, respectively. They are statistically significant for all three years (P1, P2, and P3 = 0.00 < 0.05), suggesting that assimilating satellite sea ice data effectively improves the sea ice extent/concentration forecast.
b. Sea ice thickness
Figure 2 shows the hindcast of the ensemble mean of the averaged Arctic sea ice thickness (north of 55°N, left panel) and corresponding RMSE (right panel). At the initial stage, the simulated ice thickness in CTL is too thick for all three years with respect to the observation (a combination of SMOS and CryoSat-2), and this deviation is maintained through the model integration. LESTKF-1 makes the initial ice thickness relatively thinner compared to that of CTL (orange line vs red line). After that, similar to CTL, the initial deviation maintains through the rest of the hindcast. By contrast, LESTKF-E5 improves the ice thickness simulation remarkably, in which the deviation is continuously reduced, approaching to the observation near the end of the integration (green line).
The RMSEs of CTL, LESTKF-1, and LESTKF-E5 (right panel in Fig. 2) demonstrate similar temporal evolution as the ice thickness in 2012 and 2013, although there are some small variations for RMSEs of CTL and LESTKF-1 (Figs. 2b,d). The averaged RMSEs of CTL, LESTKF-1, and LESTKF-E5 are 1.10 (1.26), 1.02 (1.09), and 0.90 (0.87) m in 2012 (2013). All p values are statistically significant (P1, P2, and P3 = 0.00 < 0.05), suggesting that the assimilation is effective to improve sea ice thickness forecasts. However, there is not much difference between LESTKF-1 and LESTKF-E5 in 2014, and the hypothesis test suggests that LESTKF-E5 is not significantly different from LESTKF-1 (P3 = 0.14 > 0.05).
As mentioned previously, the satellite-derived sea ice thickness assimilated in the CFSv2 is a combination of SMOS and CryoSat-2. The observed ice thickness averaged over the studying period shows thick ice (>1.5 m, corresponding to multiyear ice) extending from the Canadian Archipelago and north of Greenland to the central Arctic Ocean and thin ice (<1.5 m) covering large parts of the Beaufort, Chukchi, East Siberian, and Laptev Seas for all three years (first column in Fig. 3). Apparently, in 2012 and 2013, the simulated ice thickness in CTL (second column in Fig. 3) is too thick in the entire Arctic Ocean, especially for thin ice. By contrast, the pattern of the simulated ice thickness in 2014 is very different from that of 2012 and 2013 (i.e., the ice is obviously thinner in large parts of the Beaufort and Chukchi Seas and the central Arctic Ocean, which is closer to the observation). This is also reflected by the spatial distribution of RMSE with respect to the observation averaged over the studying period (not shown). In 2012 (top panel) and 2013 (middle panel), much of the Arctic Ocean is dominated by large errors, which are relatively reduced in LESTKF-1 but significantly reduced in LESTKF-E5. However, in 2014 (bottom panel), large errors only locate in the Canadian Archipelago and the East Siberian and Laptev Seas in CTL and LESTKF-1, although the error is also lowered in LESTKF-E5 but not as broad as those of 2012 and 2013. Thus, no significant difference between LESTKF-1 and LESTKF-E5 in 2014 in Fig. 2 is primarily due to the change in the distribution of thick ice simulated by the CFSv2. This suggests that the assimilation of satellite sea ice data has a strong effect when the simulated spatial pattern of sea ice thickness is different from that of the observation.
c. Evaluation with independent data
The simulated sea ice thickness from the three hindcasts is further evaluated using independent observations. First, we compare them with IceBridge airborne observations. Figure 4 shows spatial distribution of IceBridge observations and simulated sea ice thickness. Here the simulated ice thickness is interpolated on IceBridge unstructured meshes. Note that the overlapping of IceBridge samplings may occur at some grid points for each individual year, which are demonstrated in the last column in Fig. 4. The ratio between the overlapped grid points and the total is less than 0.2%, suggesting that the overlapping is negligible. In general, CTL (second column) produces much thicker ice in the Beaufort Sea and relatively thinner ice in the Canadian Archipelago and north of Greenland for all three years. This is also the case for LESTKF-1 (third column), although the bias is reduced somewhat in the Beaufort Sea. By contrast, LESTKF-E5 (fourth column) shows much better agreement with IceBridge observations in the Beaufort Sea, but there is no improvement in the Canadian Archipelago and north of Greenland.
We further analyze the temporal evolution of the daily statistics of IceBridge and simulated sea ice thickness. Here the simulated sea ice thickness for each day having IceBridge observations is interpolated from the model grid to the IceBridge track. In Fig. 5, the center of the circle represents the daily mean ice thickness and the radius of the circle represents 0.5 standard deviations of the ice thickness. In general, the simulated ice thickness underestimates the observed deviation for all three years, which is largely because the resolution of the IceBridge track is higher than the model grid. Compared to IceBridge, the simulated ice thickness with sea ice assimilation shows better agreement than that without any data assimilation, not only for the daily mean but also for the standard deviation, especially for LESTKF-E5 (green circle). However, the improvement is primarily seen before late March. After late March, the result with sea ice assimilation is even worse than that without sea ice assimilation. Figure 6 provides the difference between the simulated and IceBridge ice thickness after 26 March. During this period, IceBridge measurements are only located in the north of Greenland and/or the Canadian Archipelago where the thickest and most heavily ridged ice is present. In this area, compared to the observation, the simulated ice thickness in CTL (first column) and LESTKF-1 (second column) shows mixed positive and negative biases. By contrast, LESTKF-E5 (third column) is dominated by negative bias. This is primarily because the CryoSat-2 ice thickness assimilated into the model is thinner than IceBridge ice thicknesses as well as thinner than the CTL simulation. Nevertheless, we calculate the sea ice thickness probability density function (PDF) distribution for the simulation, IceBridge, and their difference (Fig. 7). It appears that for all three years, the PDF of LESTKF-E5 becomes tighter around IceBridge (Fig. 7, upper panel) and much more normal-looking compared to that of CTL and LESTKF-1, which deviate significantly from normality (Fig. 7, lower panel).
Second, we compare the simulated sea ice thickness with and without sea ice assimilation against a buoy measurement outside the coverage of IceBridge (Fig. 8). Again, without assimilating satellite sea ice data, the model bias is large (>1.5 m) compared to the buoy observation. LESTKF-E5 agrees best with the observation during the drifting period with a reduced RMSE of 1.03 m. This is also reflected by the PDF of the difference between the simulated and IceBridge ice thickness (not shown), suggesting that the assimilation significantly improves sea ice thickness forecast.
5. Discussion and conclusions
In this study, the ensemble-based local error subspace transform Kalman filter (LESTKF) is employed in the NCEP climate forecast model to assimilate 1) sea ice concentration from SSMIS and 2) sea ice thickness from a combination of SMOS and CryoSat-2 that provide basinwide thickness information. Three years (2012, 2013, and 2014) are selected for hindcast experiments. For each year, three numerical experiments are carried out, including 1) CTL (without sea ice assimilation), 2) LESTKF-1 (with initial sea ice assimilation), and 3) LESTKF-E5 (with every 5-day sea ice assimilation).
Our results show that compared to CTL, LESTKF-1 improves the prediction of sea ice extent/concentration and thickness initially, but the initial improvement tends to gradually diminish after ~3-week integration for the ice extent, whereas it remains quite steady through the model integration for the ice thickness. Here, we calculate the autocorrelation coefficient for the time series of the simulated total Arctic sea ice extent of CTL (without any sea ice assimilation). It decreases with increasing time, and the sign is reversed after 2–3 week integration. This suggests that the initial improvement due to sea ice assimilation may diminish in ~(2–3) weeks in the CFSv2, which might be model dependent. Walsh and Johnson (1979) examined modes of variability in Arctic sea ice area and found an autocorrelation of about 6 days up to 3 months lag for the main modes. Lukovich and Barber (2007) showed a coherent sea ice concentration persistence pattern with a time scale of 3–7 weeks.
LESTKF-E5 improves sea ice concentration and thickness forecast remarkably. The deviations in both sea ice extent and thickness in CTL are significantly reduced and approach the satellite observation by the end of the hindcast. This would benefit sea ice forecast beyond the study period. It is also noticed that assimilating satellite-derived sea ice thickness has a strong effect when the simulated spatial pattern of ice thickness is different from that of the observation. The comparison of the simulated sea ice extent and thickness for each individual 5 days is also examined since the forecast duration is 5 days for LESTKF-E5. Here the predictive skill for each individual 5 days is defined as follows:
where X = LESTKF-1 or LESTKF-E5. The closer the value to 1, the higher the predictive skill. The result shows that the predictive skill of LESTKF-E5 is higher than that of LESTKF-1 for each individual 5 days, which is particularly true for sea ice thickness (not shown).
The modeled sea ice thickness is further evaluated against independent sea ice thickness observations. Evaluation with IceBridge measurements shows that assimilating satellite sea ice data improves the ice thickness prediction in the Beaufort Sea. Evaluation with the sea ice mass balance buoy measurement also confirms significantly improved forecast of the ice thickness when satellite sea ice data are assimilated. As mentioned previously, the satellite-derived sea ice thickness assimilated in the CFSv2 is a combination of SMOS and CryoSat-2. We also carry out an additional model experiment for the year 2012. In this experiment, only SMOS sea ice thickness less than 1 m is assimilated (LESTKF-E5SMOS). In general, the evolution of the simulated sea ice extent and thickness of LESTKF-E5SMOS behaves like that of LESTKF-E5. However, compared to both the satellite and IceBridge observations, LESTKF-E5SMOS has larger bias and RMSE for sea ice thickness (as well as sea ice extent/concentration) than those of LESTKF-E5. It is noted that the assimilation of the CryoSat-2 sea ice thickness is helpful for reducing the ice thickness bias extending from the Canadian Archipelago and north of Greenland to the central Arctic Ocean simulated by the CFSv2 (not shown).
To provide additional perspective on the predictive skill, we also examine the simulated sea ice extent against the persistence prediction. The persistence prediction assumes the forecasted sea ice extent will not change through the hindcast. We calculate the ratio of sea ice extent bias (relative to the satellite observation) between three hindcasts and the persistence prediction. It appears that the ratio decreases through the model integration. The predictive skill becomes higher than that of the persistence prediction after ~(3–4) weeks of integration. This is particularly true for LESTKF-E5, which has a much smaller ratio than those of CTL and LESTKF-1 (not shown).
To further determine to what extent the aforementioned improved sea ice prediction may result from assimilating sea ice thickness from SMOS and CryoSat-2, we perform an additional model experiment for the year 2012, in which only sea ice concentration from the SSMIS is assimilated every 5 days in the CFSv2 (LESTKF-E5SIC; still nine ensemble members). As shown in Fig. 9, the simulated ice extent of LESTKF-E5SIC behaves more like that of LESTKF-1 (blue line vs orange line in Fig. 9a). The improved skill gradually diminishes after ~3-week integration, getting close to CTL, but RMSE in LESTKF-E5SIC is relatively smaller (but statistically significant) than that of LESTKF-1 (Fig. 9b). The simulated ice thickness of LESTKF-E5SIC is clearly worse than that of LESTKF-E5 (Fig. 9c), which is also evident in the averaged RMSE of sea ice thickness (Fig. 9d). Thus, the improvement of sea ice thickness prediction is minor by only assimilating sea ice concentration. By contrast, the simulation with continuous assimilation of sea ice thickness has remarkably reduced systematic errors over the hindcast compared to that with sea ice concentration assimilation only or without any data assimilation. The averaged RMSE over the studying period has been reduced from 1.08 m in LESTKF-E5SIC to 0.90 m in LESTKF-E5.
Figure 10 shows the RMSE of the simulated sea ice thickness and concentration in LESTKF-E5 and LESTKF-E5SIC with respect to the observations. It appears that large RMSEs of the simulated sea ice thickness in the Beaufort, Chukchi, East Siberian, and Laptev Seas in LESTKF-E5SIC are greatly reduced in LESTKF-E5 as the satellite ice thickness is assimilated (Fig. 10a vs Fig. 10b). Meanwhile, large RMSEs of the simulated sea ice concentration in the above areas in LESTKF-E5SIC are also reduced in LESTKF-E5 (Fig. 10d vs Fig. 10e). This suggests that assimilating satellite-derived sea ice thickness can improve the predictive skill of sea ice thickness dramatically, which also benefits the forecast of sea ice extent/concentration somehow due to the coupling between the ice thickness and concentration. For example, change in sea ice thickness modifies heat transfer from the ocean, and thinner sea ice results in increased ocean heat loss and thus further warms the atmosphere.
Recent research suggested that the persistence of sea ice thickness anomalies is much higher than that of sea ice concentration anomalies (e.g., Blanchard-Wrigglesworth et al. 2011; Wang et al. 2013; Day et al. 2014a,b). Thus, the corrected state of sea ice thickness would provide a system memory that allows for forecast with potential skill over longer periods. The improvement in sea ice prediction could be larger during the melting season when the covariance between ice thickness and concentration is strong. The SMOS and CryoSat-2 sea ice thickness therefore provide valuable information to improve thickness forecast when assimilated in sea ice model components of coupled climate forecast systems. Currently satellite-derived sea ice thickness (SMOS and CryoSat-2) is primarily available during the cold season, since the ice thickness retrieval algorithms are compromised in summer in the presence of melt. It is critical to develop datasets of year-round sea ice thickness, covering broad areas of the Arctic that would aid in the forecast procedures. A few efforts are making progress; for example, there is a new approach to estimate sea ice thickness in summer in the presence of melt pond (Istomina et al. 2016). An alternative way to deal with this is to generate sea ice thickness analysis, like the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) sea ice thickness data (Collow et al. 2015).
In the aforementioned numerical experiments, we use nine ensemble members to represent model realizations, compromising between a reasonable model spread and computational cost. To examine the sensitivity of sea ice prediction to the number of ensemble members, we conduct an additional experiment for the year 2012 using 15 ensemble members (LESTKF-E5EM15) to increase the rank of model ensemble covariance. As shown in Fig. 11, it is intuitive that a larger model spread is evident as the ensemble member increases. Over the studying period, the agreement of the forecasted ice extent of LESTKF-E5EM15 is slightly better than that of LESTKF-E5 (blue line vs green line in Fig. 11a), but the hypothesis test suggests that the RMSE of the forecasted ice concentration of LESTKF-E5EM15 is not significantly (p = 0.97 > 0.05) different from that of LESTKF-E5 (Fig. 11b). The predicted ice thickness of LESTKF-E5EM15 closely resembles that of LESTKF-E5 (Fig. 11c), but the RMSE of the predicted ice thickness of LESTKF-E5EM15 is relatively increased compared with that of LESTKF-E5 (0.94 m for the 15-member ensemble and 0.90 m for the 9-member ensemble; Fig. 11d).
Additionally, as mentioned in section 3, for the initialization step, we generate nine ensembles that provide estimates of initial model states and uncertainties prior to the evolution of the hindcast by integrating the CFSv2 for one month (e.g., from 15 February to 15 March in 2012). Here, we also generate another set of the initial ensembles by integrating the CFSv2 for three months (LESTKF-E5MON3) that includes intraseasonal sea ice variability and repeat the every-5-day assimilation experiment for the year 2012 (Fig. 12). It appears that over the studying period, the prediction of both sea ice extent (p = 0.46 > 0.05) and thickness (p = 0.65 > 0.05) of LESTKF-E5MON3 is not significantly different from that of LESTKF-E5.
In summary, our experiments demonstrate how sea ice forecast in the NCEP Climate Forecast System can be improved by assimilating satellite-derived sea ice concentration and thickness. Certainly, there are potentials for further optimization and extension—for example, using different observational errors for SMOS and CryoSat-2 or spatially varying observational errors and assimilating additional variables that influence atmosphere–sea ice–ocean interactions (e.g., sea ice velocity and sea surface temperature). These will be investigated in future studies.
We thank the NSIDC for providing SSMIS sea ice concentration (http://nsidc.org/data/nsidc-0079), the University of Hamburg for providing SMOS sea ice thickness (https://icdc.cen.uni-hamburg.de/thredds/catalog/ftpthredds/smos_sea_ice_thickness/catalog.html), meereisportal.de for providing CryoSat-2 sea ice thickness (http://data.seaiceportal.de), the NSIDC for providing the operational airborne IceBridge sea ice thickness (https://nsidc.org/data/docs/daac/icebridge/evaluation_products/sea-ice-freeboard-snowdepth-thickness-quicklook-index.html), and the ice mass balance buoy program for providing buoy-measured sea ice thickness (http://imb-crrel-dartmouth.org/imb.crrel/buoysum.htm). This research is supported by the Climate Observation and Earth System Science Divisions, Climate Program Office, NOAA, U.S. Department of Commerce (NA15OAR4310163 and NA14OAR4310216), and the NSFC (41676185).