Estimating parameters in a sea ice model using an ensemble Kalman filter

Uncertain or inaccurate parameters in sea ice models influence seasonal predictions and climate change projections in terms of both mean and trend. We explore the feasibility and benefits of applying an ensemble Kalman filter (EnKF) to estimate parameters in the Los Alamos sea ice model (CICE). Parameter estimation (PE) is applied to the highly influential dry snow grain radius and combined with state estimation in a series of perfect model observing system simulation experiments (OSSEs). Allowing the parameter to vary in space improves performance along the sea ice edge but degrades in the central Arctic compared to requiring the parameter to be uniform everywhere, suggesting that spatially varying parameters will likely improve PE performance at local scales and should be considered with caution. We compare experiments with both PE and state estimation to experiments with only the latter and have found that the benefits of PE mostly occur after the data assimilation period, when no observations are available to assimilate (i.e., the forecast period), which suggests PE's relevance for improving seasonal predictions of Arctic sea ice.

1 Introduction

Arctic sea ice has undergone rapid decline in recent decades in all seasons (e.g., Stroeve et al., 2012; Serreze and Stroeve, 2015). The frequent large deviations of Arctic sea ice cover from its climatology and the impact of sea ice cover on the overlying atmosphere and on ocean–atmosphere fluxes motivate including an active sea ice component in seasonal-to-sub-seasonal (S2S) weather forecasts (Vitart et al., 2015). The persistence and reemergence of sea ice thickness (SIT) and sea surface temperature anomalies are major sources of predictability for Arctic sea ice extent (SIE; Blanchard-Wrigglesworth et al., 2011). Previous studies have demonstrated the importance of accurate initial conditions, especially SIT, in predicting Arctic sea ice extent (Day et al., 2014). Hence studies applying data assimilation (DA) techniques to fuse observations with model simulations are actively investigated (e.g., Lisæter et al., 2003; Chen et al., 2017; Massonnet et al., 2015), most of which are focused on improving model states only, not the parameters in sea ice parameterization schemes.

Sea ice models, like other components of Earth system models, can suffer large uncertainties originating from uncertain parameters. The widely used Los Alamos sea ice model version 5 (CICE5), given its various complex schemes, has hundreds of uncertain parameters, such as in the Delta-Eddington shortwave radiation scheme (Briegleb and Light, 2007). The default values of these parameters are usually chosen based on point measurements that are taken on multi-year sea ice (Light et al., 2008). Urrego-Blanco et al. (2015) conducted an uncertainty quantification study of CICE5 and ranked the parameters based on the sensitivities of model predictions to a list of parameters. This work provides guidance on which parameters could be estimated using an objective method and during which seasons. Their findings suggest that the estimates of the Arctic sea ice area and extent are especially sensitive to certain parameters (e.g., snow conductivity and snow grain size) in summer. However, they also discussed that their sensitivities could be low as a consequence of prescribing atmospheric forcing in their model setup, so parametric uncertainties are expected to be larger year-round (particularly in winter) in a fully coupled model. Previous studies suggest that the ensemble spread of sea ice states is generally small in winter (e.g., Lisaeter et al., 2003; Fritzner et al., 2018; Zhang et al., 2018), which will lead to limited update on model state variables or parameters. Also, sea ice concentration (SIC) reaches 100 % in most of regions in winter and hence does not leave enough room for improvements by DA. The ensemble spread in summer, however, is much larger. Since we run stand-alone CICE5 given that our aim is to demonstrate the utility of parameter estimation (PE) for sea ice, we conduct DA experiments with PE in summer.

Two types of observations are assimilated in our study, sea ice concentration and thickness (SIC and SIT, respectively). Satellite-retrieved SIC observations are widely utilized in the sea ice DA community, while the application of SIT observations is more challenging given its large uncertainty and lack of data in summer (Zygmuntowska et al., 2014). Previous studies on Arctic sea ice predictability emphasized the importance of summer SIT observations (e.g., Day et al., 2014; Dirkson et al., 2017). We explore the benefits of SIT observations (in addition to SIC) on sea ice parameter estimation and advocate the need to extend the data coverage of SIT observations into late spring and summer, which is actually possible in ICESat-2 (Kwok et al., 2020).

Despite the importance of sea ice model parameters, few studies have tried to estimate or reduce the parametric uncertainties, partly due to the large effort and computational cost if parameter calibration is done in a trial-and-error fashion. A more systematic way is through DA. Anderson (2002) demonstrated the feasibility of updating parameters using an ensemble filter in a low-order model. Annan et al. (2005) were among the first to apply an ensemble filter to estimate parameters in a complex Earth system model. Massonnet et al. (2014) employed the ensemble Kalman filter (EnKF) in a sea ice model to estimate three parameters that control sea ice dynamics. In addition to achieving their goal of improving the sea ice drift, they also realized slight improvements in the SIT distribution and extent as well as in the sea ice export through the Fram Strait.

Our purpose is to expand upon previous studies to explore the feasibility of optimizing sea ice parameters by asking how different observations (concentration and thickness in this study) would constrain the parameters differently, whether we need to allow parameters to vary spatially, and what the benefits are of the updated parameters both when observations are available for assimilation (the DA period) and when observations are not available (the forecast period).

2 The sea ice data assimilation framework

We use CICE5 linked to the Data Assimilation Research Testbed (DART) (Anderson et al., 2009) within the framework of the Community Earth System Model version 2 (CESM2) (http://www.cesm.ucar.edu/models/cesm2, last access: May 2019). The ocean is modeled as a slab ocean, and the atmospheric forcing is prescribed from a DART/CAM ensemble reanalysis (Raeder et al., 2012). Details of this framework can be found in Zhang et al. (2018). The default DART/CICE framework is only used for state estimation; we extend DART/CICE to include parameter estimation in this study. During the assimilation, DART and CICE5 cycle between a DA step with DART and a 1 d forecast step with CICE5. During the DA step, the selected sea ice variables are placed into a “DART state vector” that is to be passed to the filter. The DART state vector is augmented by adding selected sea ice parameters, so that the parameters and state variables are both updated by the filter in the same way. The updated state variables are then post-processed (if needed) and sent with the updated parameters back to CICE5 for the next 1 d forecast step. The post-process step is necessary when the updated variable goes beyond its physical boundaries, for example, when SIC is negative or larger than 100 %. Unlike state variables, the parameters are not modified during CICE5 forecast steps.

3 Experiment design and evaluation methods

The parameter we select, R_snw, represents the standard deviation of dry snow grain radius that controls the optical properties of snow and is one of the key parameters that determine snow albedo in the Delta-Eddington solar radiation parameterization treatment (Briegleb and Light, 2007). We pick R_snw because it is one of the parameters that the model predictions are sensitive to (Urrego-Blanco et al., 2016) and is also one of the parameters perturbed to generate ensemble spread in Zhang et al. (2018). Instead of directly tuning snow albedo that could result in inconsistencies with the rest of the parameterization scheme, tuning R_snw changes the inherent optical properties of snow in a self-consistent fashion (Briegleb and Light, 2007). Increasing R_snw leads to smaller dry snow grain radius and larger snow albedo (Hunke et al., 2015). The default value of R_snw is 1.5, which corresponds to a fresh snow grain radius of 125 µm (Holland et al., 2012). Many parameters in CICE5, like R_snw, have default values based on limited field observations. As sea ice models increase in complexity, empirical parameters will increasingly need to be calibrated objectively. More comprehensive observations at a large scale will presumably benefit a better representation of snow and ice properties in sea ice models.

The configurations of conducted experiments are listed in Table 1. We begin with a free run of CICE5 without DA (hereafter FREE) with 30 ensemble members. Each ensemble member has a unique value of R_snw, which is constant in time and space. The ensemble of R_snw values is randomly drawn from a uniform distribution spanning −2 to 2. One of the ensemble members is designated as the truth with the true value of R_snw. Following Zhang et al. (2018), synthetic observations are created by adding random noise to SIC and SIT taken from the truth ensemble member. The noise follows a normal distribution with zero mean and a standard deviation of 15 % for SIC and 40 cm for SIT. The FREE experiment does not assimilate any observations, and the R_snw values stay the same throughout the experimental period.

Table 1 List of experiments with different configurations and RMSE of the total Arctic sea ice area and volume calculated over two experiment periods: DA (April to October 2005) and forecast (April to September 2006) for the seven experiments. All the experiments use the same localization half-width and prior inflation algorithm as stated in Sect. 3.

We then conduct two pairs of experiments to test the feasibility of estimating parameters using the ensemble adjustment Kalman filter (EAKF) (Anderson, 2002), which is a deterministic ensemble square root filter. Each experiment assimilates daily SIC or SIT synthetic observations. The first pair is referred to as DAsicPEcst and DAsitPEcst, with the former assimilating SIC observations and the latter SIT observations. In the first pair, each ensemble member has a unique spatially uniform R_snw. The second pair is referred to as DAsicPEvar and DAsitPEvar, which has a separate value of R_snw at each horizontal grid point. The augmented state has the single parameter for R_snw in the first pair or the two-dimensional grid of R_snw parameters in the second pair.

All variables in the sea ice state vector are two-dimensional in space. The parameter R_snw and the state variables were updated based on their correlations with neighboring observations. The posterior ensemble generated by DART is always spatially varying. For the first pair of experiments, we take an area-weighted average of the two-dimensional posterior to get a spatially invariant R_snw to send back to CICE5. For the second pair of experiments, the spatially varying posterior R_snw was sent to CICE5. In all experiments, the sea ice component is run for 1 d to produce a new state that is augmented with the posterior R_snw of the previous DA step. To increase the prior ensemble spread of R_snw, a spatially and temporally adaptive inflation was applied to the priors of both the model states and R_snw before they are sent to the filter (Anderson, 2007). The initial value, standard deviation, and inflation damping value of the adaptive inflation are 1.0, 0.6, and 0.9, respectively. The localization half-width is 0.01 radians (about 64 km) as discussed in Zhang et al. (2018). We also reject observations that are 3 standard deviations of the expected difference away from the ensemble mean of the forecast.

A third pair of experiments is conducted with only state DA (no parameter estimation), known as DAsic and DAsit, which assimilate daily SIC and SIT synthetic observations, respectively. DAsic and DAsit have the same ensemble set of R_snw, which is also the initial set of R_snw in the above PE experiments. The ensemble of R_snw remains fixed throughout the experiment period.

All experiments begin on 1 April 2005 and run for 18 months. Synthetic observations are assimilated only during the first 6 months (the DA period), and the next 12 months are a pure forecast period to mimic the real-world situation when making a forecast. The values of R_snw hold constant once DA ceases. We do not perform DA beyond October 2005 for two reasons. First, sea ice states have small ensemble spread in winter, as illustrated in Fig. 1a, so DA updates tend to be small. In contrast, the relatively larger spread from April to October ensures that assimilating observations can have more impact in updating model state variables and parameters. Second, the snow albedo feedback only influences the sea ice state when sunlight is present.

Figure 1 Time series of (a) the Arctic sea ice area and (b) sea ice volume from a free CICE5 run. Each gray line represents one ensemble member, black line the ensemble mean, and red line the truth. Time series of (c) the parameter R_snw for two DA experiments. Blue line represents DAsicPEcst that assimilates SIC observations, and magenta represents DAsitPEcst that assimilates SIT. The red reference line indicates the true value of R_snw. Each error bar represents 2 standard deviations of the 30 ensemble members of R_snw. Error bar is shown for every 5 d.

Several commonly used error indices are calculated to evaluate the performance of the experiments. The root-mean-square error (RMSE) of Arctic SIE and the area-weighted spatially averaged root-mean-square error (RMSE_t) are defined as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(\stackrel{\mathrm{‾}}{x_i^m}-x_i^t{)}^{\mathrm{2}}}{N}}},\\ & {\mathrm{RMSE}}_t=\sqrt{{\displaystyle \frac{{\sum }_{j=\mathrm{1}}^M(\stackrel{\mathrm{‾}}{x_j^m}-x_j^t{)}^{\mathrm{2}}}{M}}},\end{array}\end{array}$$

where i and j are the indices in time and space, x refers to Arctic SIE in RMSE and may refer to parameters or model states in RMSE_t, N is the number of days, and M is the number of grid cells. The superscripts m and t refer to model and truth, respectively. The overbar indicates the mean of the model ensemble.

Model bias is defined as the mean of the 30-member ensemble of the experiments minus the truth. Absolute bias difference (ABD) between two experiments is defined as follows:

$$\begin{array}{}\text{(2)}& \mathrm{ABD}=\left|\stackrel{\mathrm{‾}}{x_i^{\mathrm{case}\mathrm{1}}}-x_i^t\right|-\left|\stackrel{\mathrm{‾}}{x_i^{\mathrm{case}\mathrm{2}}}-x_i^t\right|,\end{array}$$

where x may refer to parameters or model states, the superscripts t refers to the truth, and case1 and case2 refer to the two experiments to compare. The overbar indicates the mean of the model ensemble.

4 Results and discussion

4.1 Temporally and spatially invariant parameters

The ensemble mean of FREE underestimates SIC throughout the year (Fig. 1a) partly because our arbitrary ensemble member selected as the truth has an above-average R_snw (Fig. 1c). As such, we would intuitively expect R_snw to have a positive increment as a result of assimilating SIC observations. Figure 1c confirms that R_snw increments are positive, with the posterior ensemble mean gradually approaching the true value during the DA period in the spatially constant PE experiments (DAsicPEcst and DAsitPEcst). The posterior R_snw has smaller ensemble spread than the prior R_snw (also see Fig. S1d, e, and f), which is consistent with the EAKF theory. In Fig. 1c DAsitPEcst outperforms DAsicPEcst starting in June, indicating that SIT provides more information than SIC for R_snw. Similarly, with state-only DA, Zhang et al. (2018) found that SIT is more efficient than SIC observations at constraining state variables. There could be several reasons why the rate at which R_snw approaches the true value decreases with time. First, the ensemble spread of R_snw may be insufficient because no uncertainty is introduced into R_snw in CICE5 during the forecast step. It is an open question how much additional uncertainty should be introduced into the parameters. To help avoid filter divergence, we apply the prior adaptive inflation to the parameters (as well as to the model states), which may still be not enough. Second, the correlation between R_snw and the observations may be too weak. Solar radiation becomes very low by the end of September, and hence R_snw has little impact on sea ice, which explains the weak correlation between R_snw and the observations (further discussed below). Either reason could result in a negligible update to R_snw.

The correlations between R_snw and the observations have unique spatial patterns and evolve with time. On 1 May, the correlation between R_snw and SIC is generally positive (Fig. 2a). The positive correlations are significant especially where SIC is under ∼100 %. Larger R_snw corresponds to higher snow albedo and more reflected sunlight, which in turn delays the melting of sea ice. The correlations are still significant along the ice edges in August (Fig. 2c) and become noisier and have less significant values by the end of the melt season (Fig. 2e). The correlation between R_snw and SIT has different spatial patterns (Fig. S2b, d, and f). Negative correlations between R_snw and SIT on 1 May can be seen in the Chukchi Sea, Beaufort Sea, and East Siberian Sea, where R_snw and SIC have positive correlations. This suggests that where SIC increases with R_snw in spring it is possible that SIT actually decreases, which might be due to elevated concentration raising the compressive strength and reducing sea ice deformation. While a brighter surface is able to reduce thickness over large regions in spring, the effect is mostly gone by the end of summer, when positive correlation prevails.

Figure 2 Correlations between (a) R_snw and SIC and (b) R_snw and SIT for 1 May 2005, (c) R_snw and SIC and (d) R_snw and SIT for 1 August 2005, and (e) R_snw and SIC and (f) R_snw and SIT for 1 October 2005. At each point, we calculate the correlation of R_snw and the observed quantities across the 30 ensemble members on the selected dates. The posterior states outputted from the experiments DAsicPEcst and DAsitPEcst are used for calculation.

4.2 Spatially varying R_snw

We discussed in Sect. 4.1 that processes relating R_snw and observed quantities have complex spatial features. The spatial map of the posterior R_snw  and the reduction in the ensemble spread of R_snw after EAKF in the first pair of experiments (Fig. S1) also suggest that the updates are concentrated on the ice marginal zones. It may be too crude to use a single value of R_snw for the whole Arctic. We let R_snw be a spatially varying parameter in the second pair of PE experiments, even though the true R_snw is spatially invariant. The spatial features of R_snw will purely depend on how R_snw correlates with the observations. As in DAsicPEcst and DAsitPEcst, the analysis field of R_snw is spatially varying, and we did a spatial averaging to get a single number for the next run. R_snw along the sea ice edges gets updated more, while R_snw in the center is less influenced. But the averaging smoothed out this spatial feature. In DAsicPEvar and DAsitPEvar, we let the spatially varying 2D analysis field of R_snw be the R_snw field in the next run, so the spatial feature was carried along the simulation.

Figure 3 depicts the ABD of R_snw (defined in Sect. 2) between different pairs of experiments at the end of the DA period. Figure 3a and d confirm that DAsicPEcst and DAsitPEcst improve the R_snw compared to FREE. Figure 3b and e show the spatial feature of improvements or degradations in R_snw for the two spatially varying PE experiments. They both show the contrast between the ice marginal zones and the central Arctic. Improvements are mostly seen along the ice edges. Spotty improvements in the inner Arctic can be found in DAsitPEvar (Fig. 3e), while degradations are prevailing in the inner Arctic in DAsicPEvar (Fig. 3b). Figure 3c and f highlight the improvements or degradations from allowing R_snw to vary spatially. The general features are that DAsicPEvar and DAsitPEvar have reduced R_snw biases more along the ice edges compared with DAsicPEcst and DAsitPEcst. However, degradations (Fig. 3c) or negligible improvements (Fig. 3f) are found in the central Arctic. This suggests that spatially invariant PE generally works better for the whole pan-Arctic regions, while spatially varying PE can work well in the ice marginal zones but not in the central Arctic, especially when SIC is the only observed quantity. SIC has little variability in the central Arctic, and hence assimilating the SIC observations will not add much information for parameters or model states. Besides the improvements along the sea ice edges, the SIT DA also has a benefit in the inner ice pack (Fig. 3e), which is consistent with the results of the first pair of experiments that SIT in general provides more information than the SIC observations, especially in the regions where SIC has little variability. However, spatially varying R_snw has small advantages over spatially invariant R_snw in the ice marginal regions but degradations in the central Arctic too (Fig. 3f). The degradations in R_snw but improvements in SIC (Fig. 5a and c; discussed in Sect. 4.3) in the central Arctic suggest that R_snw is likely over-adjusted to cancel out other errors (e.g., noise from atmospheric forcing fields).

Figure 3 The differences of absolute mean bias (ABD; see Eq. 2) of R_snw between the DA experiments – (a) DAsicPEcst, (b) DAsicPEvar, (d) DAsitPEcst, and (e) DAsitPEvar – and the control experiment FREE, and between the spatially varying PE experiments and the spatially constant PE experiments: (c) DAsicPEvar and DAsicPEcst, and (f) DAsitPEvar and DAsitPEcst.

4.3 Additional improvements in model states

We demonstrated that R_snw approaches the true value by assimilating SIC or SIT (at different rates) in the previous sections. We now investigate whether PE also improves the simulation of model states, beginning with time series of the pan-Arctic sea ice area and volume in all of our experiments (see Fig. 4).

Figure 4 Daily biases of (a) the total Arctic sea ice area and (b) the total Arctic sea ice volume for FREE (black), DAsic (blue), DAsicPEcst (green), DAsicPEvar (purple), DAsit (orange), DAsitPEcst (pink), and DAsitPEvar (red). Gray dash line in each plot represents the zero reference line. The blue line in panel (a) is overlapped by the purple and green lines in the first half of time. The black line in panel (a) is overlapped by the orange and blue lines in the second half of time. The black line in (b) is overlapped by the blue line from February to July.

In our preceding work, we showed that assimilating SIC and SIT could improve model states (Zhang et al., 2018), which can also be confirmed in Fig. 4. During the DA period, DAsic can efficiently reduce biases in area, but DAsic has limited influence on volume. Within about a month into the forecast period, DAsic improves neither area nor volume. In contrast, DAsit is highly beneficial at reducing both area and volume during the DA period, with at least some improvement to volume persisting through the whole 1-year forecast period.

We find that updating R_snw has a relatively large impact on volume beginning in spring of the forecast period (Fig. 4b). Treating R_snw as either a spatially varying or a constant parameter has about the same effect until late summer of the forecast period. In fact, all of the PE experiments outperform the state-only DA experiments in the forecast period. As shown in Table 1, SIT DA with PE always performs the best, reducing the bias in area by up to 63 % and reducing the bias in volume by up to 73 %. SIC DA with PE is second best in terms of simulating the area, reducing the bias by up to 37 %. SIC DA with PE is comparable to DAsit in simulating volume, reducing the bias by around 30 %.

Finally, we compare the spatial patterns of bias reduction in SIC and SIT from PE experiments by comparing RMSE_t of SIT in DAsicPEcst and DAsitPEcst to their state-only DA counterparts, DAsic and DAsit (see Fig. 5). The comparisons are made in two periods: the DA period (April to October 2005) and the forecast period (April to September 2006). Zhang et al. (2018) showed that the DAsic could only improve SIT along the sea ice edges. Figure 5a demonstrates that DAsicPEcst offers some improvements in the central Arctic as well. Improvements that resulted from a more accurate R_snw in the forecast period are more prominent (Fig. 5b). For DAsitPEcst, SIT is improved almost everywhere in the Arctic, with slight degradations along the ice edges (Fig. 5c). The improvements persist throughout the forecast period (Fig. 5d).

Figure 5 The relative differences of RMSE_t of SIT between DAsicPEcst and DAsic for the (a) DA experiment period and (b) forecast period, and between DAsitPEcst and DAsit for the (c) DA experiment period and (d) forecast period. The differences of RMSE_t are divided by the RMSE_t of DAsic and DAsit, respectively, to get the relative differences.

5 Conclusions

We extend the functionality of DART/CICE to do PE through the EAKF as well as updating the model states. One of the key parameters determining sea ice surface albedo, R_snw, is estimated as an example in this study. R_snw is updated using the filter. We designed a series of OSSEs to demonstrate the feasibility of PE in CICE5. Results show that R_snw gradually approaches the true value during the DA period (from April to October 2005). Updating parameters with PE could further improve the model state estimation but not prominently in the DA period. During the forecast period, with a better representation of the parameter, the PE experiments show significant superiority over the state-only DA experiments, both in SIC and SIT. The results in the forecast period indicate that, by updating parameters as well as state variables, assimilating SIC observations only is comparable to assimilating SIT observations. We concluded that SIT is the most important variable to be observed in Zhang et al. (2018), but satellite observations of SIT have large uncertainties and only cover a short time period. We could alternatively improve model parameters by assimilating SIC observations, with the ultimate goal of improving SIT. Results from the subset of experiments treating R_snw as a spatially varying parameter suggest that the R_snw biases are mostly reduced along the sea ice edges but not as much in the central Arctic. We suggest that varying R_snw spatially is not effective when conducting DA for the whole Arctic, but worth exploring when it comes to regional studies, such as in the seasonal sea ice zones.