Sea Ice modelling

This section and the next one focus on the “forward integration” spot in the centre of Figure 4.1, designing the architecture of an OOFS.

Sea ice drift forecasts are affected by multiple sources of uncertainties. The surface winds are one of the most important external forces driving the motion of the sea ice in the central Arctic (Thorndike and Colony, 1982). Moreover, the uncertainties in the atmospheric reanalysis in the Arctic are higher than those at the mid-latitudes, and observations are insufficient to estimate the statistical characteristics (scale, amplitudes) of the errors. Rabatel et al. (2018) investigated the sensitivity of sea ice drift using neXtSIM-EB for the uncertainties of the surface winds. They concluded that, in regions of highly compact ice cover, the accuracy of surface wind forcing and sea ice rheology are both important for the probabilistic forecast skill of sea ice trajectories.

Table 6.1. Overview of operating and approved satellites and sensors for the sea ice observations grouped into: ESA and Eumetsat missions (yellow), 3rd Party Missions (green) and new approved missions from ESA and NASA/CNES (blue). Spectrometers and infrared radiometers are only sensing in cloud free conditions. Note ( *) that CFOSAT flies a combined altimeter and real aperture radar at five distinct incidence angles up to 10 degrees. Harmony (**) comprises two bi-static satellites that will fly in convoy with Sentinel-1.

The ocean below the ice contains large quantities of heat and momentum, enough to melt the sea ice and to cause ice drift and deformations. Uncertainties in ocean temperature, vertical mixing, and currents are then very meaningful for the sea ice. The surface ocean salinity is important, as the melting point temperature depends on it. However, measuring ocean properties and particularly currents below the sea ice is challenging and uncertainties are rather high.

Uncertain initial conditions, particularly the sea ice thickness, persist a long time (Chevallier and Salas-Mélia, 2012). Blockley and Peterson (2018) showed that the sea ice conditions in spring persist typically a few months into the summer and are an important source of large-scale predictability. Errors in the position of the ice edge at the beginning of a forecast are usually persistent throughout the forecast run and ought to be post-processed for practical use.

Finally, sea ice models are dependent on their numerous model parameters, both in the sea ice dynamics and thermodynamics (Urrego-Blanco et al. 2016).

Initialized forecasts are critically dependent on the observations used for their initialization. To be useful for operational systems, observations are needed in near real-time for short-term forecasts and with limited time lag for seasonal and longer forecasts. There are unique challenges involved in polar observations because of its remoteness, harsh conditions, and long polar night. However, forecasting systems are making use of satellite observations for initialization, most routinely for sea ice concentration. Additionally, new products, such as sea ice thickness and drift, are becoming available and may ultimately improve the predictive capabilities. 

Sea ice reconnaissance flights were mostly occasional until after the second world war, with the exception of the USSR which started systematic flights with Polar Aviation as early as 1929 to monitor the Northern Sea Route. The USA and Japan gradually increased the frequency of their flights at the turn of the 1950’s and adopted the WMO sea ice charting standard proposed in 1952 (WMO, 1970). These flights are still used nowadays, mostly in Canada, but have elsewhere been superseded by satellite data.

Regular and routine sea ice observations are today performed by a variety of satellite-based measurements provided by several space agencies (as grouped in the matrix in Table 6.1 and organised by satellite sensor classes). This has been accomplished thanks to a large number of major technical and scientific milestones and achievements over more than 40 years, as further addressed below. Note that spectrometers and radiometers are only sensing in cloud free conditions.

Started in 1978, the longest satellite record to cover the whole Arctic comes from polar orbiting passive microwave sensors onboard the satellites SMMR, SSM/I, AMSR-E and AMSR2 (Cavalieri and Parkinson, 2012) which provide the sea ice areal concentration. Their main advantage is that they can see through clouds but still a few issues remain, especially with the summer ice, because the sensor does not properly discriminate between open water and signatures from wet snow and melt ponds. This and other technical issues are accommodated differently in a multitude of algorithms that calculate sea ice concentrations from the raw passive microwave retrievals (Ivanova et al., 2014, 2015). This is an important matter for data assimilation as we will see in Section 6.2.5. The resolution of passive microwaves depends on the frequency band used, with the most precise low-frequency channels having the largest footprint (as large as 60 km). However, gridded sea ice concentration data can be found at resolutions between 6 and 25 km, which is consistent with current operational models of the whole Arctic but still coarse with respect to the needs of any operational users navigating in ice-infested waters. SAR and satellite data in the visible channels (VIIRS, AVHRR, MODIS, SPOT) provide much more detail at spatial resolutions finer than 1 km, which is what the ship captain would need, for example to detect and sail along a lead. However, both types of data suffer from poor coverage, SAR images because the acquisition frequency may be limited, and visible data because they are impaired by the frequent cloud coverage and by Arctic winter darkness.

For short-term forecasts, it is important to assess how the sea ice is moving. Various sea ice drift products are obtained from different satellites and can be split in two types: 1) the coarse resolution, full spatial coverage products using passive microwave radiometers and scatterometers (most accurate retrievals in winter because of the aforementioned limitations of passive microwave data during summer; see review by Sumata et al., 2014); and 2) the high resolution but reduced coverage SAR-based products (Kwok 2006). The SAR coverage has recently significantly improved by the launch of the ESA Sentinel-1 A/B missions offering full daily coverage in high latitudes (Korosov and Rampal. 2017). In comparison, drifting buoys on sea ice still provide the longest (more than 40 years) data record of the IABP but with limited spatial coverage.

Table 6.2. Overview of sea ice variables observed per group of satellite sensors listed in Table 6.1. Note (*) that spectrometers and radiometers are only sensing in cloud free conditions.

Sea ice thickness observations from satellites have recently become routinely available. These use different principles to obtain either: sea ice freeboard of thick ice, for example from the satellite altimetry missions ICESat, ICESat2 (Kwok et al., 2007), CryoSat-2 (Laxon et al., 2013) and Sentinel-3; and the thickness of thin ice, derived from the SMOS (Tian-Kunze et al., 2014). These observations are quite complex and come with relatively high uncertainties (Zygmuntowska et al, 2014, Tian-Kunze et al., 2014). As discussed above, sea ice thickness is an important source of sea ice predictability on seasonal and longer timescales. Other aspects of the sea ice, such as snow cover, snow thickness and melt pond characteristics, may also be important for sea ice forecasts on seasonal and longer timescales. Remote sensing is used to characterise these aspects of the sea ice. For example, snow depth information is being provided through the NASA Operation IceBridge airborne campaign (Kurtz et al., 2013), and melt pond fractions have been derived from satellite data in the visible channels (Rösel et al., 2012). Combined use of IceSAT-2, CryoSat-2 and Altika has also demonstrated promising capabilities to recover reliable snow depth estimates during winter (Guerreiro et al., 2016). Ice mass buoys are also providing in-situ measurements of snow depth and other sea ice characteristics (Richter-Menge et al., 2006; Perovich et al., 2008). However, only limited work has been done to quantify the possible influence of these types of observations for forecasting systems.

As already indicated in Table 6.1, the continuity of sea ice observations from satellites are indeed assured by the approved future satellite missions such as CIMR, Cristal, and ROSE-L Copernicus Sentinel Expansion missions under preparation by ESA. However, more dedicated field campaigns are still needed to assess the uncertainties of the satellite-based retrievals, as well as to harvest the multi-sensor synergies as can be noted from Table 6.2. In turn, this would improve the quality and use of satellite data, and expectedly advance the forecast skill of sea ice on seasonal to inter-annual timescales.

Most modern large-scale sea ice models are based on very similar foundations. The ice is generally modelled as a continuum using a Eulerian perspective, with the sea ice moving in a horizontal plane, subject to both external and internal forces. The dynamic evolution of the sea ice cover is described using two continuity equations and the momentum equation. The thermodynamic evolution is modelled within each column of the grid and is modelled as a heat diffusion process within the slab of sea ice. There are substantially varying degrees of complexity in the treatment of the thermodynamic processes, ranging from treating all the ice as being of a single thickness (Hibler, 1979) and treating the heat diffusion without resolving the temperature profile (Semtner, 1976), to using multiple thickness categories (Hibler, 1979, and numerous later variations) and treating the heat diffusion using mushy-layer dynamics (Feltham, et al., 2006). 

The main equations for a simple dynamic model of sea ice with two categories (ice and open water) are the two continuity equations and the momentum equation. The continuity equation for mass is:

Where m is the sea ice mass per unit area, S_m is a thermodynamic source/sink term and v is velocity. In the case of a single sea ice category the continuity equation for the sea ice distribution takes the basic form:

With A is the sea ice concentration and S_A is a source/sink term. In addition, the condition A≤1 is imposed. This can be interpreted as a ridging condition since m can increase even if A does not (Hibler, 1979). Together these equations describe the advection of the sea ice in a given velocity field.

The momentum equation is generally written as (Connolley et al., 2004):

Here k̂   is a unit vector normal to the surface, τ_a and τ_w are the air andwater stresses, f is the Coriolis parameter, g is the gravitational acceleration, ∇H is the gradient of the sea surface height and σ is the sea ice stress tensor. The acceleration term on the left-hand side may be set to zero, depending on the model implementation. The last term on the right hand side ∇·σ, describes forces due to internal stress while the other terms are all external factors. Wind and water stresses are generally treated as quadratic drag (e.g., McPhee, 1975). In the absence of internal stress, the sea ice is in “free drift” and the model simplifies drastically. Free drift forecasts have therefore been used for a long time (Grumbine 1998) are still used operationally.

The thermodynamic equation is the heat diffusion equation:

Where ρc is the heat capacity of sea ice or snow and k is the heat conductivity. This equation can be solved in various ways (see Figure 6.3) (e.g., Maykut and Untersteiner, 1971; Semtner, 1976; Bitz and Libscomb, 1999; Winton, 2000; Huwald et al., 2005), discretized in the vertical (Figure 6.3). These take into account different physical properties and numerical solutions in solving the equation.

In addition to these grid-scale quantities, many models consider various sub-grid scale information and parameterisations. The most important of those is arguably the ice thickness distribution (ITD). This assumes that each grid cell of the model contains not is of uniform thickness, but of varying thicknesses described by an ice thickness distribution g (see Figure 6.3). This is in principle a continuous distribution of thicknesses, which is modified through dynamic and thermodynamic processes. The governing equation of evolution of the ice thickness distribution is (e.g., Thorndike et al., 1975):

Where f  is the thermodynamic growth or melt rate, h is the ice thickness, and 𝚿 is a mechanical redistribution function.

In practice, sea ice models must use a discretized version of the ice thickness distribution, resulting in models with several distinct thickness categories (Figure 6.4 right, for a top view of a grid cell). The thickness distribution then becomes (Bitz et al., 2001):

Where M is the number of thickness categories, M_i is the thickness of category i, and δ(h) is the Dirac delta function. Various implementations exist, but the one from Bitz et al., (2001) with five thickness categories remains a popular choice.

In addition to these two basic components, a large number of sub-grid scale processes can and should be represented, depending on the use cases for each model. These include simulation of melt points (Flocco et al., 2010; Hunke et al., 2013), changes in atmospheric and oceanic drag due to sea ice roughness (Tsamados et al., 2014), and salt rejection from freezing sea ice (Vancoppenolle et al. 2009).

The relationship between the internal stress and resulting deformation is referred to as rheology. Basically, all continuum, geophysical-scale sea ice models currently employ the VP rheology proposed by Hibler (1979) or some direct descendant of that work. The VP rheology treats the sea ice as a continuum and assumes it deforms in a viscous manner with a high viscosity until the internal stress reaches a plastic threshold, determined by a yield curve which usually has an elliptic shape (see Figure 6.5). Several important improvements have been made to the original VP rheology (e.g. Hunke and Dukowicz, 1997; Lemieux et al., 2010; Bouillon et al., 2013; Kimmritz et al., 2017), but the physical principles remain the same.

The VP rheology has enjoyed tremendous success and is used for time scales from days to centuries and spatial scales from tens of kilometres to basin scales. However, the VP rheology is not without faults when it comes to both the underlying assumptions (see in particular Coon et al., 2007) and the results it produces. In model inter-comparison studies, there is generally a very large spread - well beyond observed variability - in key prognostic variables such as sea ice thickness, concentration, and drift (Chevallier et al., 2017; Tandon et al., 2018). The sharp gradients in velocities, which are known as LKFs that are related to ridge and lead formation, are also poorly reproduced in any VP-based model running at a coarser resolution than about 2 km - a resolution that is an order of magnitude higher than the observational data (Spreen et al., 2017; Hutter et al., 2019).

Therefore, several authors, such as Tremblay and Mysak (1997), Wilchinsky and Feltham (2004), Schreyer et al. (2006), Girard-Ardhuin and Ezraty (2012), Dansereau et al. (2016), and Ólason et al. (2022), have suggested alternative approaches to the VP rheology. The EAP rheology of Wilchinsky and Feltham (2004) was implemented in the CICE model (Tsamados et al., 2013) and has been used in several studies, although it was not widely adopted yet (Bouchat et al., 2021; Hutter et al. 2021). The brittle rheologies of Girard et al. (2011), Dansereau et al. (2016), and Ólason et al. (2022) have all been implemented in the neXtSIM model (Bouillon and Rampal 2015; Rampal et al., 2019; Ólason et al., 2022) and used for forecasting and scientific research by the team involved in the model. The current neXtSIM version uses the BBM rheology of Ólason et al. (2022).

Practically all sea ice models used in modern forecasting platforms are based on the principles described above. They use a Eulerian reference frame and use some version of the VP or the Elastic-Viscous-Plastic (EVP) rheologies. The thermodynamic growth and melt of ice are described through the diffusion of heat between ocean and atmosphere, through the ice. As such, they all follow the same general design philosophy. The main differences exist in the form of different choices of parameterisation and differences in data assimilation approaches.

Today, the CICE model (e.g. Hunke et al., 2021) is likely the most widely used sea ice model for operational forecasts. This model was developed at the Los Alamos National Laboratory and was originally designed to be part of the CCSM. Thanks to its clean and modular design, the model has been used in other multiple modelling systems, as a stand-alone model, part of sea ice-ocean models, and part of climate and earth-system models. The LIM and SI3 models (Rousset et al., 2015), which are part of the NEMO modelling system, are also very widely used, but only within the NEMO modelling system. Other sea-ice models include the SIS (Adcroft et al., 2019), which is part of the GFDL ocean modelling system MOM, the MITgcm sea-ice model (Losch et al., 2010) and the FESOM sea-ice model (FESIM, Danilov et al., 2015).

In contrast, only a few stand-alone sea ice models have used moving Lagrangian coordinates (Hopkins 2004), among which the neXtSIM model (Rampal et al., 2016; Rampal et al., 2019; Ólason et al., 2022) and the DEMSI model (Turner et al., 2022). neXtSIM-F is unique among forecasting models as it uses both a moving Lagrangian mesh and a newly developed brittle rheology, the Brittle Bingham-Maxwell (Williams et al., 2021; Ólason et al., 2022). This setup gives results that are clearly different from the classical systems, and arguably more realistic (Rampal et al., 2016; Rampal et al., 2019; Ólason et al., 2020; Ólason et al., 2022). The key improvement is a much more realistic representation of the deformation statistics of the ice cover, which gives more realistic leads and ridges in the model. Sea ice drift simulated by neXtSIM is also very realistic, and the pan-Arctic ice-thickness distribution is also quite good (Williams et al., 2021).

Sea ice models are integral parts of Earth system models. The reason for this is that at high latitudes sea ice insulates the relatively warm ocean from the cold atmosphere. Over an unbroken sea ice cover, the atmosphere can therefore cool much more than it could if it was not insulated by the presence of sea ice. This has an impact on all ocean-atmosphere interactions in the polar regions, and therefore a global climate or Earth system model without a sea ice model can simply not function.

Sea ice interacts with the atmosphere through heat, moisture, and momentum exchanges. In summer incoming shortwave radiation melts the ice surface but would warm up the ocean surface in the absence of sea ice. In winter, heat conduction from the ocean and through the ice only results in a very modest amount of heat flux to the atmosphere. However, the dominant heat flux source is radiant cooling through long wave radiation from the surface. This happens because surface cooling through long wave radiation is much more efficient than the heat conduction through ice from the ocean, resulting in a surface that is colder than the lowest layers of the atmosphere. The result is a predominant temperature inversion and a stable atmospheric boundary layer. This reduces even further the latent and sensible heat fluxes from the surface. However, openings in the ice (leads and polynyas) expose the relatively warm ocean surface to the atmospheric boundary layer, which causes mixing and breaks down the stable boundary layer. 

Momentum transfer between ice and atmosphere happens through wind stress at the surface of the ice. This is the main driver of ice movement and exerts a drag on the atmosphere, slowing down the wind. The amount of momentum transferred between ice and atmosphere is determined primarily by the stability of the atmospheric boundary layer (Gryanik and Lüpkes, 2017) and the roughness of the ice. While parameterisations and studies on the ice surface roughness have been proposed (Lüpkes et al., 2012, Castellani et al., 2014), consistent and basin-scale observations of the atmospheric drag coefficient over sea ice are currently unavailable (Petty et al., 2017). In a modelling context, our ability to predict ice surface roughness is severely limited, as most ice-atmosphere coupled models do not take surface roughness into account when calculating atmosphere-ice momentum exchanges.

Sea ice interacts with the ocean through heat, fresh-water, and salt exchanges, as well as momentum exchanges. During summer, the mixed layer may warm up due to shortwave heating through openings in the ice. This makes the ice melting from below, causing release of both freshwater and salt into the ocean. In winter, the atmosphere extracts heat from the ocean through the ice, causing new ice to form at the bottom of the existing ice pack. This causes a net heat and freshwater flux out of the ocean. However, most of the salt present in the ocean cannot enter the ice, since the ice is much fresher than the ocean (ca. 15 vs 30 PSS for newly formed ice in the Arctic). This results in a layer of very salty water forming below the ice, which then sinks into the mixed layer. The resulting salt plumes generally reach the bottom of the halocline but may also be mixed into the mixed layer in the presence of turbulence.

Momentum transfer between ice and ocean happens through interface stress between ocean and ice. The momentum coupling of ice and ocean is much stronger than that of ice and atmosphere, and the ice can be considered as the first layer in the ocean’s Ekman spiral. Ice-ocean stress drives most geostrophic flows in ice covered areas, as well as some larger scale circulation patterns.

It is also worth mentioning the interaction between sea ice and ocean waves. Waves entering the ice pack may mechanically fracture it into smaller sea ice floes. This can widen the MIZ, which may also be viewed as the area where the ice is fractured by waves. Fracturing the ice into smaller floes increases the mobility of the ice cover, the momentum transfer between atmosphere, ocean, and ice, and this may cause enhanced melting of the ice through lateral melt. The ice in turn dampens the waves causing an attenuation of the wave amplitude, so that waves will only penetrate a limited distance into the ice pack, depending on the size of the waves and the compactness of the pack. Wave-ice interactions are of major importance in the Southern Ocean, but less so in the Arctic, where much less of the ice edge is exposed to open ocean.

Virtually all climate or earth-system models today include sea-ice models of the classical description above, i.e., a Eulerian reference frame, VP or EVP rheology, and thermodynamics and column physics of varying complexity. They generally include very simplistic formulations for the momentum transfer between atmosphere, ocean, and ice, and no ice-wave interactions. This is true for all the CMIP6 models. In fact, the sea ice models used in today's forecasting models were all designed for climate modelling, the only current exception is the above mentioned neXtSIM model. It is not clear how this lineage of the models affects the quality of their short-term predictions. It could be argued that a good large-scale sea ice model should be able to represent scales from ca. 1 km up to the basin scales and from hours to centuries. This is not the current case, but the discussion of why it is this way and how to address it is still in its infancy (Hunke et al., 2020; Blockley et al., 2020).

In nearly all forecasting platforms, the sea ice model is coupled to an ocean model. There are platforms that use fully coupled atmosphere-ocean-sea ice models and only a few platforms using a stand-alone sea ice model. The reasons for this are partly historical: most sea ice models are written as parts of ocean models. Ocean forecasting and re-analysis platforms have tended to include a sea ice model from the start, making a dedicated sea ice forecasting platform redundant. In addition, the coupling between sea ice and ocean is quite strong, so running a separate sea ice forecasting platform can bring its own set of challenges. On the other hand, a stand-alone sea ice forecasting platform can be run at a higher resolution and can be used as a technology preview, as in the case of the neXtSIM-F platform.

Probabilistic forecasts, which are widely used in weather forecasting (Molteni et al., 1996; Leutbecher and Palmer, 2008), are still in their infancy in sea ice forecasting. Probabilistic predictions rely on an ensemble of model simulations (e.g. a Monte Carlo simulation) used to describe the forecast uncertainty stemming from errors in the model parameters, initial and boundary conditions, and any external forcing. The resulting range of model outputs is used to retrieve statistical information, such as the ensemble mean and its spread (i.e. the standard deviation), which are thus used instead of the deterministic forecast to estimate the associated uncertainty (see Figure 6.6). The multiple simultaneous sources of errors usually make the forecast accuracy of the ensemble mean exceed that of the single deterministic prediction (Leith, 1974), although the spread often underestimates the actual forecast error when the sources of error are not all adequately accounted for (Buizza et al., 2005). Monte Carlo techniques are already common practice in different areas (e.g. Dobney et al., 2000; Hackett et al., 2006; Breivik and Allen, 2008; Melsom et al., 2012; Motra et al., 2016; Duraisamy and Iaccarino, 2017) and a common tool for sensitivity analysis.

As introduced in the previous section, a sea ice forecast needs to regularly assimilate operational observations, which, at present, are mostly satellite data. The most tempting way forward is to insert directly the satellite observed concentrations and thicknesses into the model. However, this is not as easy as it sounds in a complex sea ice code where a large number of model variables are dependent on each other. Hence, various data assimilation methods are used for sea ice models, similar to those used for ocean physical, biogeochemical models or weather models. The most common method is nudging, which is less disruptive than direct insertion: the data are introduced gradually over a given time scale (Lindsay and Zhang, 2006). The nudged model is then assumed to adjust itself progressively using its own equations. But how much can we rely on such adjustments?

When the ocean mixed layer is too warm to sustain sea ice but observations show the presence of sea ice, a data assimilation system updating only sea ice would add sea ice on top of the warm waters, but the huge heat capacity of the ocean would then melt the added sea ice almost immediately. The ocean mixed layer temperature and salinity must be adjusted accordingly. This suggests that when used in a coupled ice-ocean system, assimilation of sea ice observations ought to be coupled in the sense that it should update both the sea ice and the ocean properties consistently. In data assimilation jargon, this means that the sea ice observation should be projected down to the ocean column using a multivariate forecast error covariance matrix.

Dynamical model ensembles are a practical way to estimate the covariances mentioned above. In data assimilation terminology, the state vector must include all prognostic variables of the coupled model (ocean and sea ice variables) and the ensemble of model runs can be used to calculate empirically the cross-covariances between sea ice and ocean variables. Similarly, observations of the ocean are used to update sea ice variables, although this situation is less common. Using an EnKF (see section 5.5.2), Lisæter et al. (2003) demonstrated that the coupled assimilation of sea ice properties can modify the ocean surface temperatures in rather systematic ways (adding sea ice cools down the water), but not ocean salinities. However, according to sea ice halodynamics, the freezing of sea ice injects salty brines to the ocean mixed layer and the melting releases fresher water, but these simple relationships explain only a part of the sea ice-salinity cross-covariances and a relationship may arise in other situations without the intervention of sea ice thermodynamics: the wind may occasionally blow the sea ice on top of more saline water. Sakov et al. (2012) showed how the sea ice-salinity cross-covariance can change sign on either side of the ice edge in the Barents Sea: the sea ice-salinity correlation turns negative on the ocean side because the main process responsible for melting is the advection of warm and saline Atlantic water near the surface, thus the sea ice-salinity correlation is made through the intermediate of the surface temperature variable. The last finding does not hold in locations where the sea ice is isolated from the Atlantic water, but such isolation may not remain forever if the open water mixing reaches these warm waters (Rippeth et al., 2015). The assimilation of sea ice concentrations with the EnKF described in Lisæter et al. (2003) was included in the near-real-time TOPAZ forecasts in 2003.

An alternative to ensemble methods is the use of an adjoint model as in the 4D-variational (4D-Var) data assimilation method. The adjoint model and the tangent linear model calculate the sensitivity of observed variables to the control variables within the duration of the assimilation window. If tangent linear and adjoint models are available both for the ocean and the sea ice models, they can exchange information about the interface variables, like the heat, salt, and momentum fluxes. Since these correlations are usually monovariate at the beginning of the assimilation window, the length of the assimilation window should be as long as possible. The most recent experiments report successful applications of the 4D4D-VarVAR in an Arctic regional configuration for durations of one year or longer (Fenty and Heimbach 2013; Fenty et al., 2015); an adjoint model for the EVP sea ice rheology has been introduced later (Toyoda et al. 2019). The advantage of the 4D-Var method is that it returns one optimised model trajectory, which is very useful for oceanographic interpretation (Kauker et al., 2009) and quantitative network design (Kaminski et al., 2015) but, to our knowledge, 4D-Var is not used for operational sea ice-ocean forecasting.

A computationally simplified variant of 4D-Var is known as 3D-Var, in which the same increment is used to compute the model equivalent of the observation-minus-reference state differences at all times in the assimilation. Owing to the relatively low cost of the scheme compared with the full 4D-Var, 3D-Var is commonly used by operational forecasting centres around the world (Usui et al., 2006; Mogensen et al., 2012; Hebert et al., 2015; Tonani et al., 2015; Waters et al., 2015; Smith et al., 2016, see Table 6.3 below).

Coupled multivariate covariances do not necessarily cure all the troubles of assimilating sea ice observations. Another source of problems is the lack of respect of the traditional Gauss-linear assumptions underlying classical data assimilation methods. By definition, sea ice concentrations have bounded values between zero and one, while other sea ice tracer variables (thickness, snow depth) have positive values only. Ocean temperatures are not allowed below the freezing point. While it should be easy for a monovariate assimilation method, based on a heuristic covariance function, to preserve monotonicity and therefore the bounds of variables (Wackernagel, 2003), an ensemble-based covariance (or a tangent linear model) may generate values out of bounds. Honouring the bounds can be forced by different means, either by nonlinear transformations of the variables (a method called Gaussian anamorphosis in geostatistics; Bertino et al., 2003, Barth et al., 2015) or by including inequality constraints in the cost function (Lauvernet et al., 2009; Simon et al., 2012; Janjic et al., 2014). Altogether, the benefits of multivariate flow-dependent covariances still outbeat the inconvenience of values out of bounds.

There are continuous improvements to data assimilation methods in chaotic high-dimensional systems, such as coupled sea ice-ocean models. But new models and new observations always call for further developments in data assimilation. In particular, sea ice models expressed in Lagrangian grids with automatic remeshing are uncommon targets for data assimilation. Ensemble Kalman Filtering techniques rely on cross-covariances between observed and unobserved variables, which implies that the grid cells have to be uniquely identified across different members of the ensemble. This also becomes difficult when adaptive remeshing is turned on, unless the Lagrangian model output is interpolated on a fixed grid, at the risk of smoothing the very localised kinematic features (long cracks, ridges and leads) that they are meant to simulate (Aydoǧdu et al., 2019). Lagrangian models do not offer any easy differentiation/automatic adjoint capabilities, thus preventing the use of variational techniques. It should also be noted that a coupling framework such as CESM is sufficiently flexible to allow several instances of a model component to be run (for example, the atmosphere) for each instance of another (for example, the ocean), thus allowing to use different data assimilation methods for the sea ice, ocean, land, and atmosphere. An important aspect in view of coupled data assimilation and ensemble forecasting is that the uncertainties are consistent across these compartments; the error statistics at the base of the atmosphere are consistent with those at the surface of the sea ice and similarly between the bottom of the sea ice and the ocean surface. This is possible to enforce if all components of the coupled system use an ensemble to represent the errors.

Since a measure of RMS errors of sea ice concentrations depend on arbitrary choices made by the person doing the scoring (these errors diminish as more open ocean is included in the validation area), more targeted sea ice validation metrics express the skill as distance of the forecast from the observed ice edge (Dukhovskoy et al., 2015). In the Arctic, the skill of the 24-hour forecast of ice edge location is about 50 km for the TOPAZ (including both seasonal biases and RMS errors, updated at http://cmems.met.no/ARC-MFC/) and 40 km down to 30 km depending on the input data sources in the ACNFS (Hebert et al. 2015; Posey et al., 2015), although both methods may differ and be sensitive to special configurations of the ice edge. The area of discrepancy is accepted as an objective metric with the IIEE introduced by Goessling et al. (2016). The dependence of such metrics on spatial scales can be further included in the evaluation using the FSS (Melsom et al., 2019) and an extension of the IIEE metric proposed by Goessling and Jung (2018) for the evaluation of ensemble forecasts of the ice edge. Examples of IIEE and FSS are shown in Figure 6.7.

It should be noted that the metrics related to isolines (like the ice edge, classically defined at the 15% ice concentration isoline, or at other critical values such as 50% and 85%) apply also to other isolines, like the frontier between FYI and MYI, or to theMIZ/pack boundary. Contingency tables are also a valuable approach to the validation of sea ice concentrations (Smith et al., 2016), as well as the threat scores or the Heidke Skill Score.

The forecast skills for sea ice drift have received comparatively less attention but errors in sea ice drift are important, both for their contribution to the displacement of the sea ice edge and for their cumulative contribution to the sea ice thickness distribution. Long climate simulations hint for a seasonal dependence of the forecast skills, also noted in free drift simulations (Grumbine 1998). Biases of sea ice drift have been revealed in IPCC simulations (Tandon et al., 2018) related to the seasonal cycle. To our knowledge, there are no signs that these shortcomings are corrected in recent forecast models (Xie et al., 2017, for the TOPAZ system; Hebert et al., 2015 for the ACNFS), although a review of global reanalysis systems shows that some models simulate correctly the minimum sea ice drift in March (Chevallier et al., 2017). Hebert et al. (2015) also noted that the forecast of drifter positions beats persistence although the forecast of drift speed does not, indicating that the drift direction is better forecasted than the drift speed. How to remedy these shortcomings? Although adjusting the mean speed of sea ice can be easily achieved by tuning the drag coefficients, there is no simple tuning that can make the sea ice accelerate over years or shift its seasonal cycle.

The assimilation of sea ice drift data has been so far less successful than that of sea ice concentrations: Stark et al. (2008) showed a 50% reduction of errors in ice speed but no benefit to ice concentrations and Sakov et al. (2012) indicated a low sensitivity of the sea ice drift to external perturbations in the wind forcing, which points to a shortcoming of the TOPAZ4 version of the EVP sea ice rheology. Qualitatively, the large-scale patterns of sea ice drift can be reproduced by such a model (see a typical situation in Figure 6.8) but the observed gradients between areas of low sea ice drift (North of Greenland) and strong sea ice drift (North of the Barents Sea) are smoothed by the model, which tends to simulate intermediate values of the sea ice drift speed. The forecast of 24-hours ice trajectories and locations exhibits an RMS error of 6.3 km in TOPAZ4 (Melsom et al., 2015), which does not seem to beat a simple free drift predictor (5 km in Grumbine, 1998). It should be noted that the validation is done against different data sources (sea ice drift from satellite SAR images versus IABP buoys) and at different periods (years 2012-2015 versus 80’s and 90’s decades). The SIDFEx (Goessling et al., 2020), carried out in the framework of the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) ice camp, has been the first to collect forecasts from international systems and has shown that a consensus forecast could be successfully used to order detailed satellite images of the ice camp in advance. Beyond the use of RMS errors, several alternative metrics for sea ice drift validation have been reviewed by Grumbine (2013). Validation metrics for an ensemble of trajectories from a probabilistic ice drift forecast have been proposed by Rabatel et al. (2018) and refined in Cheng et al. (2020) based on an analogy with Search and Rescue operations, in which the ensemble of trajectories define a search ellipse; the success of the forecast is the probability of containment of the target inside the ellipse.

The forecast of sea ice thickness also suffers from excessive smoothness: thick sea ice is too thin and thin sea ice is too thick (Johnson et al., 2012) and errors reach easily one to two metres. There is a dynamical contribution to these errors with the too high sea ice drift speed North of Greenland exaggerating the transport of thick sea ice into the Beaufort Gyre. However, thermodynamic contributions cannot be excluded either (in particular from snow and melt ponds). More generally, any error in the model initial conditions, atmospheric and ocean boundary conditions or its inherent parameters will eventually accumulate in sea ice thickness biases, which means that different errors can cancel each other and yield a correct sea ice thickness for the wrong reasons. It is worth stressing the important role of snow on sea ice as an effective insulator, its presence can inhibit both the growth and melt of sea ice and thus reduce its seasonal cycle. Snow predictions in sea ice-ocean models are very dependent on the quality of precipitation from weather analyses and forecasts which are difficult to validate and usually vary from one product to another (Lindsay et al., 2014).

Information on formats and types of outputs by all kinds of OOFS can be found in Chapter 4. In this Section, we list the variables related to sea ice forecasts:

• sea ice concentration (SIC)
• sea ice thickness (SIT)
• sea ice drift velocity in x- and y- directions (SIUV)
• snow depth (SNOW)
• sea ice age
• sea ice albedo (SIALB)
• sea ice temperature 

Sea ice forecasting systems generally comply with CF standards. The CF metadata conventions are a widely used standard for atmospheric, ocean, and climate data. Standard names are defined in a CF Standard Name Table (see 🔗4 ). Standard variable names from the CMIP nomenclature can be found in Notz et al. (2016).

Most present day short-term forecast systems (listed in Table 6.3) assimilate sea ice concentration and are therefore expected to perform well at forecasting the ice edge. These systems include the Canadian RIOPS (Smith et al., 2021), the United States ACNFS/GOFS3.1 (Hebert et al., 2015), the Italian GOFS16 (Iovino et al,. 2016) the Global and the Arctic Marine Forecasting System (TOPAZ, Sakov et al. 2012) by the Copernicus Marine Services. Stand-alone sea ice models, like neXtSIM-F (Williams et al., 2021), are also used for forecasting purposes and, given that their control vector excludes the ocean, they can be initialised more flexibly than coupled ice-ocean systems. Baltic forecasting systems are omitted for brevity. Ocean data assimilated are also omitted from Table 6.3.

Table 6.3.  List of present-day short-term global and Arctic forecast systems. Note that the output spatial resolution may differ from the native resolution.

* Note that the resolution of a Lagrangian triangular mesh is not comparable to square grids, thus the output resolution is 3 km.  
** VENUS is deployed on demand.