Assimilating summer sea-ice concentration into a coupled ice–ocean model using a LSEIK filter

Abstract

The decrease in summer sea-ice extent in the Arctic Ocean opens shipping routes and creates potential for many marine operations. For these activities accurate predictions of sea-ice conditions are required to maintain marine safety. In an attempt at Arctic sea-ice prediction, the summer of 2010 is selected to implement an Arctic sea-ice data assimilation (DA) study. The DA system is based on a regional Arctic configuration of the Massachusetts Institute of Technology general circulation model (MITgcm) and a local singular evolutive interpolated Kalman (LSEIK) filter to assimilate Special Sensor Microwave Imager/Sounder (SSMIS) sea-ice concentration operational products from the US National Snow and Ice Data Center (NSIDC). Based on comparisons with both the assimilated NSIDC SSMIS concentration and concentration data from the Ocean and Sea Ice Satellite Application Facility, the forecasted sea-ice edge and concentration improve upon simulations without data assimilation. By the nature of the assimilation algorithm with multivariate covariance between ice concentration and thickness, sea-ice thickness fields are also updated, and the evaluation with in situ observation shows some improvement compared to the forecast without data assimilation.

1. Introduction

In the past 30 years Arctic sea-ice extent and volume have consistently decreased in all seasons, but the maximum decline is observed in summer. In the past 10 years this summer decline has accelerated, with record lows in September. The rate of decrease in September ice extent recorded from 1979 to 1998 was (0.032 ±0.017) x 10⁶km² a^–1, and from 1999 to 2010 it increased further to (0.154±0.038) x 1 0 ⁶ km² a^–1 (Cavalieri and Parkinson, 2012; Stroeve and others, 2012). From 2003 to 2008, the observed basin-wide decline of Arctic sea-ice thickness reached 0.17 m a^–1 (Kwok and others, 2009). The melt season extended by almost 20 days from 1979 to 2007 (Markus and others, 2009). According to the latest climate model predictions, the September ice extent will drop to 1.7 x 10⁶km² in the mid-2040s and an ice-free state will be reached in 2054-58 in high-emission scenarios (Liu and others, 2013). In such a rapidly changing Arctic, the Arctic shipping routes are expected to be open in the near future and Arctic maritime activities will become more and more frequent. These activities urgently require accurate sea-ice real-time forecasts (Eicken, 2013).

In support of the CHInese National Arctic Research Expeditions (CHINARE), a simple Arctic sea-ice-ocean forecasting system based on the Massachusetts Institute of Technology general circulation model (MITgcm) (Marshall and others, 1997; see Section 2) was designed at the National Marine Environmental Forecasting Center of China (NMEFC; Yang and others, 2012). Sea-ice concentration data were incorporated with a simple reinitialization scheme. In this scheme, initial sea-ice concentration was replaced by satellite observations, and sea-ice thickness and sea surface temperature were adjusted accordingly. Although the scheme makes full use of the concentration data, it is dynamically crude and introduces physical inconsistencies. For example, Yang and others (2011), after applying the scheme, observed reinitialization shocks that lead to unrealistic sea-ice extents. While this system behavior may also be related to systematic deficiencies in the model configuration or atmospheric forcing, it is plausible to assume that forecasts and system behavior can be improved by replacing the simple reinitialization by more sophisticated data assimilation techniques that allow the combination of different types of observations and the model in a smooth, systematic way.

Several studies have demonstrated the feasibility and the benefit of assimilating observed sea-ice concentration into coupled ice–ocean models. Lisæter and others (2003) used an ensemble Kalman filter (EnKF) to assimilate Special Sensor Microwave/Imager (SSMI) concentration. Lindsay and Zhang (2006) employed a nudging scheme to assimilate monthly averaged concentration. Stark and others (2008) used an optimal interpolation method to assimilate the SSMI concentration. Wang and others (2013) developed a combined optimal interpolation and nudging scheme to assimilate concentration. In all of these studies, the assimilation of observed ice concentration in ice–ocean models was shown to improve the simulated concentration, but the improvement in ice thickness was always small. In a recent study, Tietsche and others (2013) assimilated ice concentration observations with a simple Newtonian relaxation into a coupled climate model and updated the sea-ice thickness using a proportional dependence between concentration and mean thickness. This simple scheme with the implicit assumption of a correlation between ice thickness and concentration was found to be successful in correcting sea-ice concentration and thickness. However, the model physics were crudely parameterized by a homogeneous local correlation in their analysis.

In this study, a local ensemble-based singular evolutive interpolated Kalman (SEIK) filter (Pham and others, 1998; Pham, 2001) was used to assimilate the sea-ice concentration into MITgcm over 3 months in summer. This fully dynamic SEIK filter includes the full correlations between ice thickness and concentration based on ensemble model simulations. The effectiveness of the assimilation is analyzed by comparing to the assimilated ice concentration data and a different satellite observation product. In addition, the influence of the assimilation on the ice thickness is assessed with in situ measurements.

2. Model and Atmospheric Forcing

The sea-ice module within the MITgcm (Marshall and others, 1997) includes state-of-the-art dynamics and simple zero-layer thermodynamics (Losch and others, 2010). It has been used in regional Arctic studies at varying resolution (Losch and others, 2010; Nguyen and others, 2011, 2012). We use a regional MITgcm configuration similar to those in our Arctic modeling and forecasting experiments. The modeling domain covers a limited Arctic area with open boundaries near 558 N in the Atlantic and Pacific sectors. A global configuration (Menemenlis and others, 2008) is used to provide monthly boundary conditions for potential temperature, salinity, and current velocities. The grid covering the Arctic domain is locally orthogonal and has a variable horizontal resolution with an average spacing of 18 km. The sea-ice and ocean equations are solved on the same horizontal mesh. The vertical resolution is greatest in the upper ocean, with 28 vertical levels in the top 1000 m. Bathymetry is derived from the US National Geophysical Data Center (NGDC) 2 min global relief dataset (ETOPO2) (Smith and Sandwell, 1997). The monthly mean river runoff is based on the Arctic Runoff Data Base (ARDB; Nguyen and others, 2011). The dynamics of the MITgcm sea-ice model is based on a variant of the viscous–plastic (VP) dynamic– thermodynamic sea-ice model (Zhang and Hibler, 1997) with momentum equations solved implicitly on a C-grid (Losch and others, 2010). The coupled sea-ice–ocean model is stepped forward synchronously with a time step of 1200 s.

We illustrate the development of the assimilation scheme for our Arctic model system with the help of historical forecasts (i.e. ‘hindcasts’) with analysis data (Japan Meteorological Agency Climate Data Assimilation System (JCDAS)). These analysis data start in January 2005. They are consistent with the data assimilation used in the Japanese 25 year reanalysis (JRA-25) (Onogi and others, 2007).

3. Data Assimilation Approach and Observation Data

The simulated sea-ice concentration and satellite-derived sea-ice concentration are combined using a sequential SEIK filter with second-order exact sampling (Pham, 2001) as coded within the Parallel Data Assimilation Framework (PDAF; Nerger and Hiller, 2013; http://pdaf.awi.de). The SEIK filter algorithm has been demonstrated to have some advantages over the other filters: for example, it is better suited for nonlinear models, computationally more efficient, and more accurate than the EnKF (Nerger and others, 2005). The SEIK filter has already been successfully used to assimilate sea-ice motion in a stand-alone sea-ice model (Rollenhagen and others, 2009). The filter algorithm can be divided into four phases: initialization, forecast, analysis and reinitialization.

The data assimilation process is initialized by an optimized ocean-sea-ice spin-up run (Nguyen and others, 2011). The initial uncertainties in the model ice concentration and thickness are approximated by an ensemble generated from a multivariate empirical orthogonal function (EOF) analysis of the model dynamics under variable atmospheric forcing (see Losa and others, 2012, 2014). For simplicity, the initial state error covariance matrix of the sea-ice concentration and sea-ice thickness is estimated based on a model integration over the period 1 June to 31 August 2010. In a real application, one would use a similar sampling period from the previous model year, or maybe even averaged over many previous model years. One time slice per day is collected into a set of 92 state vectors. Each state vector includes maps of sea-ice thickness and concentration. Together these 92 state vectors form a matrix that is decomposed into EOFs. The leading EOFs are used to generate an ensemble of initial ice concentration and thickness and, therefore, to approximate the background forecast error covariance. After this initialization phase, the ensemble evolves dynamically in time, driven by atmospheric forcing, to produce a forecast at the time when new data are available (here 24 hours). In the analysis and reinitialization step, the ensemble forecast is combined with the observations to create updates of the background error covariance and the ensemble states based on model-data misfit and the error statistics. After the analysis step, the reinitialized ensemble members are again propagated by the model to produce the next ensemble forecast. For more details the reader is referred to Nerger and others (2006).

Following Nerger and others (2006) and Janjic and others (2011) the SEIK analysis is applied locally at each model gridpoint given the model forecast and observational information only within a certain radius. In our study, a radius of 7 gridpoints (-126 km) was introduced to localize the analysis. To account for missing data around the North Pole (the ‘polar gap’ due to the inclination of the remote-sensing satellite), the localization radius was gradually increased from 8 gridpoints (-144 km) at 86° N to 29 gridpoints (-232 km) at the North Pole. This approach allows us to extrapolate the observed information from the surrounding regions. Within the localization radius, the observations are weighted according to their distance from the center gridpoint (Hunt and others, 2007) by a fifth-order polynomial function that mimics a Gaussian distribution (Gaspari and Cohn, 1999). Due to the statistical nature of the analysis update and the incomplete sampling of the error covariance, the analysis step may generate overshoots of too small (negative) and too high values of ice concentration and thickness. Negative ice concentrations and thicknesses are locally replaced by zero, and ice concentrations are bounded from above by 1 (100% ice cover). An ice-concentration-dependent adjustment of the ice thickness is also applied at the same time. Zero ice thickness in the presence of nonzero ice concentration is set to a new ice thickness of 2 m times ice concentration (Tietsche and others, 2013). Any ice thickness with zero concentration is removed. The sea surface temperature is not part of the state vector and hence is not modified directly by the SEIK filter; in the presence of sea ice, sea surface temperatures are updated implicitly by the model assumption of thermodynamic equilibrium between sea ice and the ocean surface water layer.

Two types of daily sea-ice concentration data are used in this study. The ice concentration observations used in the assimilation are derived from US Defense Meteorological Satellite (DMSP) F-17 Special Sensor Microwave Imager/ Sounder (SSMIS) passive microwave data, processed by the NSIDC with the NASA Team algorithm (Cavalieri and others, 2012; http://nsidc.org/data/docs/daac/nsidc0051_gsfc_ seaice.gd.html). These data are interpolated to the model grid. On average, an estimate of uncertainty in the observed sea-ice concentration is -10% (Tonboe and Nielsen, 2010), but since the errors of satellite-derived sea-ice concentration are far larger in summer than in winter (Comiso and others, 1997), and to account for a representative error, a constant value of 0.30 is used for the uncertainties in NSIDC SSMIS sea-ice concentration for summer 2010. The ice concentration data used for evaluation are from the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) Ocean and Sea Ice Satellite Application Facility (OSISAF) (Eastwood and others, 2011). The final product consists of daily fields provided on a 10km polar stereographic grid. Note that OSISAF concentration for summer 2010 is derived from another passive microwave sensor, SSM/I on board DMSP F-15, so it is independent of the NSIDC data used in the assimilation.

Satellite-based observations of ice thickness are a challenge (Kwok and Sulsky, 2010), and at present there are very few reliable summer sea-ice thickness products available. Instead of remote-sensing data we compare our simulation and assimilation results to measurements of sea-ice draft from the Beaufort Gyre Experiment Program (BGEP) upward-looking sonar (ULS) moorings located in the Beaufort Sea (BGEP_2009A, BGEP_2009D; Fig. 1) and sea-ice thickness data obtained from autonomous ice mass-balance buoys (IMBs; Perovich and others, 2009) as independent observation data for ice thickness. The error in ULS measurements of ice draft is estimated as 0.1 m (Melling and others, 1995). Drafts are converted to thickness by multiplying by a factor of 1.1, which is approximately the ratio of mean sea-water density of 1024 kg m^–3 to sea-ice density of 910 kg m^–3 (Nguyen and others, 2011). Two acoustic rangefinders on the IMBs monitor the position of the ice bottom and the snow/ice surface. The sea-ice thickness is estimated from these positions. The accuracy of both sounders is 5 mm (Richter-Menge and others, 2006). In this study, the IMB_2010A and IMB_2010B were used; their trajectories during summer 2010 are shown in Figure 1.

Fig. 1. Locations of sea-ice thickness observation and buoy trajectories from 1 June to 3 1 August 2010: BGEP_2009A (magenta square), BGEP_2009D (red square), IMB_2010A (blue line) and IMB_2010B (green line).

The performance of any data assimilation system depends on the prior model and data error statistics and how these statistics evolve in time (Losa and others, 2012, 2014). In the application here, the most crucial parameters of the data assimilation and forecasting system are the prior observational data errors, the localization length scale and inflation of the time-evolved forecast error statistics. A series of sensitivity tests has been carried out to calibrate our DA system. Here the results of the calibrated system implementation as described above are presented.

In this study, the system’s forecasting skills are evaluated with a series of 24 hour forecasts in which the local SEIK (LSEIK) filter is applied every day at 00.00 UTC over the period 1 June to 31 August 2010. The summer (June to August) of 2010 was chosen as our experimental period to check the performance of the assimilation system. The atmospheric circulation between June and August 2010 was characterized by the Arctic dipole anomaly, an atmospheric pressure pattern that contributed to the record sea-ice loss in 2007 (Wang and others, 2009). The summer of 2010 was the first time open water was found in the interior pack ice near the North Pole as early as 12 July (NSIDC, http://nsidc.org/arcticseaicenews/2010/07/.

4. Results

At each analysis step, sea-ice thickness and concentration are updated based on the available data and the forecast error covariance. The accumulation, i.e. the sum, of these update increments from 2 June to 31 August 2010 is shown in Figure 2. The spatial distribution of the accumulated increments shows that there is a systematic tendency to overestimate the ice concentration and thickness in the coastal seas that the filter algorithm tries to correct.

Fig. 2. Accumulated analysis increments of (a) sea-ice concentration and (b) sea-ice thickness (m) over the period 2 June to 3 1 August 2010. The increments refer to the update during the analysis.

Figure 3 shows the effect of assimilating NSIDC SSMIS concentration data on the simulated sea-ice concentration for 7 June (Fig. 3a and b) and 31 August 2010 (Fig. 3c and d). The strong overestimation of sea-ice concentrations in the model without data assimilation (Fig. 3a and c) is corrected towards observed values, especially in the marginal ice zone.

Fig. 3. The forecast skill improvement of sea-ice concentration on 7 June (a, b) and 3 1 August 2010 (c, d). MITgcm only (a, c) and LSEIK 24 hour forecast (b, d) minus NSIDC SSMIS ice concentration.

Figure 4 compares the temporal evolution of the root-mean-square error (RMSE) of the ice concentration forecast with and without data assimilation with respect to the assimilated NSIDC SSMIS data (Fig. 4a) and the independent OSISAF concentration (Fig. 4b) for 1 June to 31 August 2010. Following Lisæter and others (2003), all RMSEs are evaluated only at gridpoints where either the model or the observations have ice concentrations larger than 0.05. The green curve represents the RMSEs without assimilation, while the blue curve and the dots are those obtained with the LSEIK filter applied every 24 hours. As expected, the effect of the data assimilation reduces the deviations of the forecasted ice concentration from the assimilated SSMIS concentration but also from the independent OSISAF concentration. For reference, the rms difference between the OSISAF and SSMIS concentration product is shown as the black curve in Figure 4b. Although both products are derived from the similar passive microwave sensors, the deviations almost reach the RMSE of the assimilated model simulation. The 91 day forecast based on LSEIK analysis on 1 June (magenta line in Fig. 4) illustrates that updated sea-ice concentration and thickness allow for an improved forecast over a long period (much more than 5 days). For 2 months the RMSE of this forecast is smaller than the model simulation without data assimilation.

Fig. 4. Temporal evolution of RMSE differences between NSIDC SSMIS (a) and OSIAF ice concentration data (b) and MITgcm forecast (green), 9 1 day forecast based on LSEIK analysis on 1 June (magenta), mean of 24 hour ensemble forecast based on LSEIK analysis (blue), and LSEIK analysis (red) over the period 1 June to 3 1 August 2010. The deviation between NSIDC SSMIS and OSIAF concentration data is also shown as black line in (b). Date format is dd/mm.

The comparison of ice thickness predictions with in situ ULS observations (BGEP_2009A, Fig. 5a; BGEP_2009D, Fig. 5b) suggests an improvement in the sea-ice thickness with assimilation of ice concentration. Note that the numerical model carries mean thickness (volume over area) as a variable. The mean thickness is divided by the local ice concentration to arrive at the thickness shown in Figure 5. Both forecasts with and without data assimilation reproduce the gradual decrease of ice thickness. Without DA, the thickness flattens out, and the further decrease after late July is not properly simulated. In general, the agreement between predicted and observed sea-ice thicknesses has been improved by assimilating ice concentrations: the rms difference between the forecast and observations has been reduced from 0.90m to 0.57m at BGEP_2009A and from 0.95 m to 0.53 m at BGEP_2009D.

Fig. 5. Sea-ice thickness evolution at (a) BGEP_2009A, (b) BGEP_2009D Beaufort Sea, (c) IMB_2010A and (d) IMB_2010B from 1 June to 31 August, 2010: observation (black), MITgcm forecast without DA (green curve), and mean of ensemble forecast based on 24 hourly analysis (blue). Date format is dd/mm.

The drifting ice-mass balance buoys IMB_2010A and IMB_2010B are also used for comparison. The modeled sea-ice thickness is interpolated to each of the time-evolving IBM locations. For IMB_2010A (Fig. 5c), the model starts with a large positive bias of ∽ 0 . 9m on 1 June, but both forecasts capture the observed downward trend. In August, however, the assimilated model thickness increases erroneously, so the RMSE of the 3 month thickness increases from 0.83 m to 1.09m with data assimilation. The increase is caused by the analysis updates of the LSEIK filter (not shown) and can be understood as follows: In August IMB_2010A passes through the polar gap around the North Pole where no ice concentration observations are available. Instead of local data, the filter extrapolates information of low ice concentration from the fringe of the polar gap into the region and the update of thickness is purely based on prior correlations between concentration and thickness obtained from model simulations (see Section 3). The model mean thickness divided by (low) concentration as shown in Figure 5 then overestimates the observed in situ thickness. Clearly, the large data gap around the North Pole poses a limitation to our DA system, and forecasts in this region are not reliable.

Since its snow sounder failed on 7 May, the ice thickness at IMB_2010B (Fig. 5d) had to be computed from ice profile data that were available only once a week, so there are only ten data points in the period 6 June to 8 August. Similar to IMB_2010A, both model forecasts have a positive bias of ∼1.3 m on 6 June, but both forecasts capture the decreasing trend since 11 July. The LSEIK forecast is closer to the observations over most of the period, so thickness RMSE improved from 1.01 m without DA to 0.79 m with DA.

5. Conclusions

A LSEIK filter has been applied to assimilate observed sea-ice concentration data into a regional Arctic ice–ocean model. For the three summer months June to August 2010, the agreement of the ice concentration forecast with satellite observations improved in comparison with the regular model run without DA. The corrections in the mean state of the sea-ice concentration and thickness lead to an improved concentration forecast over a long period. The summer sea-ice thickness was also improved, most likely due to the multivariate covariance between ice concentration and thickness that was used in the LSEIK filter. There are, however, limits to the quality of the forecasts. The polar gap in the concentration data renders the forecasts near the North Pole unreliable and leads to obvious problems in the thickness forecasts. Still, given the fact that observed ice thickness fields are not available over the entire Arctic area in summer, our study is a step towards improving future operational sea-ice thickness forecasts.

In this study, we have improved Arctic summer sea-ice forecasts by implementing a LSEIK filter with dynamic error evolution. However, there are many other factors that can affect the forecasting behavior. Here only the sea-ice concentration observation data were assimilated, and only the ice concentration and thickness belonged to the control vector in the assimilation. The ultimate goal of a comprehensive data assimilation system would involve a full multivariate assimilation, in which variables such as ice thickness, ice-drift velocity and sea surface temperature are also updated during the analysis step of the filter algorithm. Furthermore, the model tends to overestimate the multi-year sea-ice thickness in the central Arctic (Fig. 5c and d), while the agreement with ULS data in the western Beaufort Sea is much better (Fig. 5a and b). Thickness biases were also observed by Nguyen and others (2011) who were able to reduce these biases by adjusting model parameters. Although our simulations start from their configuration, the remaining thickness biases make it clear that there is still room for improvement. Our assimilation experiments already provide some of this improvement, but it is foreseeable that including internal model parameters into the state vector and adding more data, in particular thickness data as they become available, will lead to even better agreement with observations. It remains to be seen to what extent such an experiment can contribute to better understanding of the internal model physics, their deficits, and to the improvement of model parameterizations.