# Multivariate State And Parameter Estimation With Data Assimilation Applied To Sea-Ice Models Using A Maxwell Elasto-Brittle Rheology

Yumeng Chen, Polly Smith, Alberto Carrassi, Ivo Pasmans, Laurent Bertino, Marc Bocquet, Tobias Sebastian Finn, Pierre Rampal, Véronique Dansereau

► To cite this version:

Yumeng Chen, Polly Smith, Alberto Carrassi, Ivo Pasmans, Laurent Bertino, et al. Multivariate State And Parameter Estimation With Data Assimilation Applied To Sea-Ice Models Using A Maxwell Elasto-Brittle Rheology. The Cryosphere, 2024, 18, pp.2381-2406. 10.5194/tc-18-2381-2024. insu-04604398

HAL Id: insu-04604398

https://insu.hal.science/insu-04604398

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multivariate state and parameter estimation with data assimilation applied to sea-ice models using a Maxwell elasto-brittle rheology

Yumeng Chen<fnref n="1" />, Polly Smith<fnref n="1" />,†, Alberto Carrassi<fnref n="1" />,<fnref n="2" />, Ivo Pasmans<fnref n="1" />, Laurent Bertino<fnref n="3" />, Marc Bocquet<fnref n="4" />, Tobias Sebastian Finn<fnref n="4" />, Pierre Rampal<fnref n="5" />, and Véronique Dansereau<fnref n="6" />

**Correspondence:** Yumeng Chen (yumeng.chen@reading.ac.uk)

Received: 8 August 2023 – Discussion started: 16 October 2023

Revised: 26 February 2024 – Accepted: 21 March 2024 – Published: 14 May 2024

**Abstract.** In this study, we investigate the fully multivariate state and parameter estimation through idealised simulations of a dynamics-only model that uses the novel Maxwell elasto-brittle (MEB) sea-ice rheology and in which we estimate not only the sea-ice concentration, thickness and velocity, but also its level of damage, internal stress and cohesion. Specifically, we estimate the air drag coefficient and the so-called damage parameter of the MEB model. Mimicking the realistic observation network with different combinations of observations, we demonstrate that various issues can potentially arise in a complex sea-ice model, especially in instances for which the external forcing dominates the model forecast error growth. Even though further investigation will be needed using an operational (a coupled dynamics-thermodynamics) sea-ice model, we show that, with the current observation network, it is possible to improve both the observed and the unobserved model state forecast and parameter accuracy.

## 1 Introduction

An accurate representation of the state of sea ice in the Arctic is important for making both short-term–seasonal and long-term climate predictions. Recent observations show that its extent is in decline and, in particular, that it is shifting from a multi-year ice type to “younger” first-year ice (Meier, 2017). This ongoing shift induces a larger year-to-year and inter-annual variability in the sea-ice extent, which makes the short-term and seasonal Arctic sea-ice forecasting even more challenging, despite its crucial relevance for shipping routes and fisheries (Mioduszewski et al., 2019). The Arctic sea ice is also a major player of the climate systems via its feedback to the Earth surface albedo, ocean and atmosphere global circulations. Accurate sea-ice simulations are therefore essential for better climate projections (Bertino and Holland, 2017).

As in other areas of environmental prediction, errors in sea-ice models can be attributed to errors either in the initial conditions or in the model representation of physical processes. Due to the low degree of internal instabilities of the sea-ice evolution on long space scales and timescales (> seasonal timescales), model error is the main source of prediction uncertainty, and it is especially detrimental to long-term forecasts needed for climate simulations. Numeri-

<footnote n="1">¹Department of Meteorology and National Centre for Earth Observation, University of Reading, Reading, RG6 6ET, UK</footnote>

<footnote n="2">²Department of Physics and Astronomy “Augusto Righi”, University of Bologna. Bologna, Italy</footnote>

<footnote n="3">³Nansen Environmental and Remote Sensing Center, 5007 Bergen, Norway</footnote>

<footnote n="4">⁴CEREA, École des Ponts and EDF R&D, Île-de-France, France</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, 38000 Grenoble, France</footnote>

<footnote n="6">⁶Institut des Sciences de la Terre, Institut de Géophysique de l'Environnement, Université Grenoble Alpes, CNRS, Grenoble, France</footnote>

<footnote>†deceased, 8 July 2023</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Multivariate state and parameter estimation with data assimilation applied to sea-ice models using a Maxwell elasto-brittle rheology

Yumeng Chen<fnref n="1" />, Polly Smith<fnref n="1" />,†, Alberto Carrassi<fnref n="1" />,<fnref n="2" />, Ivo Pasmans<fnref n="1" />, Laurent Bertino<fnref n="3" />, Marc Bocquet<fnref n="4" />, Tobias Sebastian Finn<fnref n="4" />, Pierre Rampal<fnref n="5" />, and Véronique Dansereau<fnref n="6" />

**Correspondence:** Yumeng Chen (yumeng.chen@reading.ac.uk)

Received: 8 August 2023 – Discussion started: 16 October 2023

Revised: 26 February 2024 – Accepted: 21 March 2024 – Published: 14 May 2024

**Abstract.** In this study, we investigate the fully multivariate state and parameter estimation through idealised simulations of a dynamics-only model that uses the novel Maxwell elasto-brittle (MEB) sea-ice rheology and in which we estimate not only the sea-ice concentration, thickness and velocity, but also its level of damage, internal stress and cohesion. Specifically, we estimate the air drag coefficient and the so-called damage parameter of the MEB model. Mimicking the realistic observation network with different combinations of observations, we demonstrate that various issues can potentially arise in a complex sea-ice model, especially in instances for which the external forcing dominates the model forecast error growth. Even though further investigation will be needed using an operational (a coupled dynamics-thermodynamics) sea-ice model, we show that, with the current observation network, it is possible to improve both the observed and the unobserved model state forecast and parameter accuracy.

## 1 Introduction

An accurate representation of the state of sea ice in the Arctic is important for making both short-term–seasonal and long-term climate predictions. Recent observations show that its extent is in decline and, in particular, that it is shifting from a multi-year ice type to “younger” first-year ice (Meier, 2017). This ongoing shift induces a larger year-to-year and inter-annual variability in the sea-ice extent, which makes the short-term and seasonal Arctic sea-ice forecasting even more challenging, despite its crucial relevance for shipping routes and fisheries (Mioduszewski et al., 2019). The Arctic sea ice is also a major player of the climate systems via its feedback to the Earth surface albedo, ocean and atmosphere global circulations. Accurate sea-ice simulations are therefore essential for better climate projections (Bertino and Holland, 2017).

As in other areas of environmental prediction, errors in sea-ice models can be attributed to errors either in the initial conditions or in the model representation of physical processes. Due to the low degree of internal instabilities of the sea-ice evolution on long space scales and timescales (> seasonal timescales), model error is the main source of prediction uncertainty, and it is especially detrimental to long-term forecasts needed for climate simulations. Numeri-

<footnote n="1">¹Department of Meteorology and National Centre for Earth Observation, University of Reading, Reading, RG6 6ET, UK</footnote>

<footnote n="2">²Department of Physics and Astronomy “Augusto Righi”, University of Bologna. Bologna, Italy</footnote>

<footnote n="3">³Nansen Environmental and Remote Sensing Center, 5007 Bergen, Norway</footnote>

<footnote n="4">⁴CEREA, École des Ponts and EDF R&D, Île-de-France, France</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, 38000 Grenoble, France</footnote>

<footnote n="6">⁶Institut des Sciences de la Terre, Institut de Géophysique de l'Environnement,
Université Grenoble Alpes, CNRS, Grenoble, France</footnote>

<footnote>†deceased, 8 July 2023</footnote>
cal models incorporate parametrisations to represent physical
processes that are not explicitly described or resolved in the
model. The parameters involved can be estimated in different
ways. One way is to base the parameter on values obtained
from laboratory experiments in idealised environments, but
this has the limitation that the selected values may not scale
up for complex realistic simulations. Alternatively, one can
select a value from a range of possible candidates based on
which candidate produces the best fit between the model and
observations (Dansereau, 2011; Miller et al., 2006; Heorton
et al., 2019). Nevertheless, as the number of parameters to
be tuned increases, the latter approach becomes computa-
tionally infeasible. In that case, an optimised version of the
aforementioned process called “data assimilation” must be
pursued.

Data assimilation (DA) combines observations with the
model forecast to provide the most likely estimate for the true
state and/or parameter of the system (see e.g. Carrassi et al.,
2018; Evensen et al., 2022; Park and Zupanski, 2022). It has
been proved to be essential for sea-ice forecasting (Sakov
et al., 2012a; Lea et al., 2015; Zuo et al., 2019). Data assim-
ilation can account for both the forecast and the observation
errors in the state and/or parameter estimation. To this ex-
tent, error statistics need to be specified. In ensemble DA,
the forecast error statistics are obtained from an ensemble of
model runs (Evensen, 2003).

Previous studies have shown that assimilating sea-ice
concentration (SIC) can reduce errors in sea-ice thickness
(SIT) (Massonnet et al., 2015) and that, in coupled sea-ice
and ocean models, sea-ice observations help in initialising
ocean fields such as the ocean and wind forcing (Toyoda
et al., 2016). Model parameters can be made part of the
state, and when it is assumed that they are affected by er-
rors, DA can be used to estimate them. An example of such
an approach in sea-ice models can be found in Massonnet
et al. (2014), who used the ensemble Kalman filter (EnKF;
Evensen, 2003) to estimate the air drag coefficients, ocean
drag coefficients and sea-ice strength parameter by assimi-
lating sea-ice drift data.

Sea-ice models include multiple variables. However, only
three of these are usually observed and assimilated: SIC,
SIT and sea-ice drift. These observations are mostly obtained
from satellites. In situ observations are usually only used for
evaluation purposes (Jakobson et al., 2012). The satellite ob-
servations of SIC show good spatial and temporal coverage
with relatively low uncertainty. Assimilation of SIC shows
tremendous benefits for sea-ice forecasting (Lisæter et al.,
2003; Stark et al., 2007). In the past decades, though less
accurate than SIC, SIT has started to be assimilated, with
further improvements to forecasts (Xie et al., 2018; Blockley
and Peterson, 2018; Fiedler et al., 2022b). The assimilation
of sea-ice drift has been less successful and has motivated
further developments of sea-ice rheological models (Sakov
et al., 2012a). These studies together show the important role
of DA in sea-ice prediction.

In this study, we explore the capability of ensemble DA to
estimate both the state and the key sea-ice parameters in a
model endowed with a Maxwell elasto-brittle (MEB) rheo-
logy. This rheology is adopted by the neXt-generation Sea Ice
Model (neXtSIM), which runs operationally on a Lagrangian
grid that uses dynamical remeshing (Rampal et al., 2016).
Due to the Lagrangian framework, the advection process
does not require explicit calculations, and the highly multi-
scale sea-ice features, especially localised features such as
the so-called linear kinematic features (LKFs), which con-
centrate sea-ice fracturing and high deformations rates, can
be easily preserved (Bouillon and Rampal, 2015). One down-
side of this framework, however, is that, with the changing
number of model grid points in time, it poses challenges
for ensemble DA. Consequently, DA approaches have been
specifically designed (e.g. Aydoğdu et al., 2019; Sampson
et al., 2021) and implemented in the context of neXtSIM
by Cheng et al. (2023). Another drawback of using a La-
grangian framework in a sea-ice model is that coupling to
existing ocean and atmospheric model components, which
are virtually all Eulerian, requires the use of a coupler on a
fixed exchange grid (Boutin et al., 2023).

Because such a coupling is essential for climate simu-
lations, a new version of neXtSIM based on an Eulerian
framework is being developed under the Scale-Aware Sea Ice
Project (SASIP; https://sasip-climate.github.io/, last access:
7 May 2024). It is in the light of these developments that
this study is inscribed. In particular, the dynamics-only sea-
ice MEB model of Dansereau et al. (2016, 2017) has been
designed following an Eulerian approach and a discontin-
uous Galerkin treatment of advection. With the aid of this
dynamics-only sea-ice MEB model, we study in detail the
capabilities, limitations and adaptations of ensemble DA to
infer state and model parameters based on synthetic observa-
tions of the Arctic sea ice.

We intentionally use a dynamics-only model whereby
thermodynamics processes are missing. Indeed, such a model
is already complex enough and, on short (daily and weekly)
timescales, sufficient to focus on how the mechanical/dy-
namical processes lead to the emergence of complex, poten-
tially nonlinear, relations between model states, model pa-
rameters and the observable quantities. Those are the rela-
tions the DA has to rely on. Thus, although the model does
not capture all of the processes at play in sea ice, we shall
see how our experiments reveal a number of complex interac-
tions. On the other hand, the use of a simpler idealised model
allows us to conduct a fully multivariate estimate where we
infer all model fields with only a handful of available ob-
servations, a realistic situation that has not yet been studied
extensively in sea-ice DA.

The paper is organised as follows. In Sect. 2 we introduce
the iterative ensemble Kalman filter (IEnKF), state-of-the-
art DA approach used in this study. Then we describe the
model and its configurations in Sect. 3. The ensemble DA
setup and twin experiments are given in Sect. 4. Results for
cal models incorporate parametrisations to represent physical
processes that are not explicitly described or resolved in the
model. The parameters involved can be estimated in different
ways. One way is to base the parameter on values obtained
from laboratory experiments in idealised environments, but
this has the limitation that the selected values may not scale
up for complex realistic simulations. Alternatively, one can
select a value from a range of possible candidates based on
which candidate produces the best fit between the model and
observations (Dansereau, 2011; Miller et al., 2006; Heorton
et al., 2019). Nevertheless, as the number of parameters to
be tuned increases, the latter approach becomes computa-
tionally infeasible. In that case, an optimised version of the
aforementioned process called “data assimilation” must be
pursued.

Data assimilation (DA) combines observations with the
model forecast to provide the most likely estimate for the true
state and/or parameter of the system (see e.g. Carrassi et al.,
2018; Evensen et al., 2022; Park and Zupanski, 2022). It has
been proved to be essential for sea-ice forecasting (Sakov
et al., 2012a; Lea et al., 2015; Zuo et al., 2019). Data assim-
ilation can account for both the forecast and the observation
errors in the state and/or parameter estimation. To this ex-
tent, error statistics need to be specified. In ensemble DA,
the forecast error statistics are obtained from an ensemble of
model runs (Evensen, 2003).

Previous studies have shown that assimilating sea-ice
concentration (SIC) can reduce errors in sea-ice thickness
(SIT) (Massonnet et al., 2015) and that, in coupled sea-ice
and ocean models, sea-ice observations help in initialising
ocean fields such as the ocean and wind forcing (Toyoda
et al., 2016). Model parameters can be made part of the
state, and when it is assumed that they are affected by er-
rors, DA can be used to estimate them. An example of such
an approach in sea-ice models can be found in Massonnet
et al. (2014), who used the ensemble Kalman filter (EnKF;
Evensen, 2003) to estimate the air drag coefficients, ocean
drag coefficients and sea-ice strength parameter by assimi-
lating sea-ice drift data.

Sea-ice models include multiple variables. However, only
three of these are usually observed and assimilated: SIC,
SIT and sea-ice drift. These observations are mostly obtained
from satellites. In situ observations are usually only used for
evaluation purposes (Jakobson et al., 2012). The satellite ob-
servations of SIC show good spatial and temporal coverage
with relatively low uncertainty. Assimilation of SIC shows
tremendous benefits for sea-ice forecasting (Lisæter et al.,
2003; Stark et al., 2007). In the past decades, though less
accurate than SIC, SIT has started to be assimilated, with
further improvements to forecasts (Xie et al., 2018; Blockley
and Peterson, 2018; Fiedler et al., 2022b). The assimilation
of sea-ice drift has been less successful and has motivated
further developments of sea-ice rheological models (Sakov
et al., 2012a). These studies together show the important role
of DA in sea-ice prediction.

In this study, we explore the capability of ensemble DA to
estimate both the state and the key sea-ice parameters in a
model endowed with a Maxwell elasto-brittle (MEB) rheo-
logy. This rheology is adopted by the neXt-generation Sea Ice
Model (neXtSIM), which runs operationally on a Lagrangian
grid that uses dynamical remeshing (Rampal et al., 2016).
Due to the Lagrangian framework, the advection process
does not require explicit calculations, and the highly multi-
scale sea-ice features, especially localised features such as
the so-called linear kinematic features (LKFs), which con-
centrate sea-ice fracturing and high deformations rates, can
be easily preserved (Bouillon and Rampal, 2015). One down-
side of this framework, however, is that, with the changing
number of model grid points in time, it poses challenges
for ensemble DA. Consequently, DA approaches have been
specifically designed (e.g. Aydoğdu et al., 2019; Sampson
et al., 2021) and implemented in the context of neXtSIM
by Cheng et al. (2023). Another drawback of using a La-
grangian framework in a sea-ice model is that coupling to
existing ocean and atmospheric model components, which
are virtually all Eulerian, requires the use of a coupler on a
fixed exchange grid (Boutin et al., 2023).

Because such a coupling is essential for climate simu-
lations, a new version of neXtSIM based on an Eulerian
framework is being developed under the Scale-Aware Sea Ice
Project (SASIP; https://sasip-climate.github.io/, last access:
7 May 2024). It is in the light of these developments that
this study is inscribed. In particular, the dynamics-only sea-
ice MEB model of Dansereau et al. (2016, 2017) has been
designed following an Eulerian approach and a discontin-
uous Galerkin treatment of advection. With the aid of this
dynamics-only sea-ice MEB model, we study in detail the
capabilities, limitations and adaptations of ensemble DA to
infer state and model parameters based on synthetic observa-
tions of the Arctic sea ice.

We intentionally use a dynamics-only model whereby
thermodynamics processes are missing. Indeed, such a model
is already complex enough and, on short (daily and weekly)
timescales, sufficient to focus on how the mechanical/dy-
namical processes lead to the emergence of complex, poten-
tially nonlinear, relations between model states, model pa-
rameters and the observable quantities. Those are the rela-
tions the DA has to rely on. Thus, although the model does
not capture all of the processes at play in sea ice, we shall
see how our experiments reveal a number of complex interac-
tions. On the other hand, the use of a simpler idealised model
allows us to conduct a fully multivariate estimate where we
infer all model fields with only a handful of available ob-
servations, a realistic situation that has not yet been studied
extensively in sea-ice DA.

The paper is organised as follows. In Sect. 2 we introduce
the iterative ensemble Kalman filter (IEnKF), state-of-the-
art DA approach used in this study. Then we describe the
model and its configurations in Sect. 3. The ensemble DA
setup and twin experiments are given in Sect. 4. Results for
both multivariate state estimation in a “perfect-model” setup and parameter estimation in a biased model are reported in Sect. 5. In Sect. 6 we discuss some of the critical aspects of ensemble DA in this context together with how to address them in future works. We finally summarise our findings in Sect. 7.

## 2 Data assimilation: the iterative ensemble Kalman filter

The iterative ensemble Kalman filter (IEnKF; Sakov et al., 2012b; Bocquet and Sakov, 2012, 2014) is a variant of the Kalman filter (KF). Like the KF, it is constructed based on the Bayes theorem under the assumption that the prior (forecast), evidence (observation) and posterior (analysis) follow a Gaussian distribution. It has successfully been applied to joint state and model parameter estimation problems (Bocquet and Sakov, 2013; Haussaire and Bocquet, 2016; Bocquet et al., 2021). In ensemble DA methods, the linear model assumption of the KF is relaxed and the prior distribution is estimated from a finite ensemble of model forecasts. A relevant feature of the IEnKF is that it solves for the analysis via a nonlinear optimisation aimed at maximising the a posteriori probability distribution. This key aspect makes it worth investigating its performance in the context of predicting Arctic sea ice, which is characterised by strong and complex nonlinear relations as well as weak nonlinearity in the observations. The IEnKF is the filter version of the more general iterative ensemble Kalman smoother (Bocquet and Sakov, 2014), and it is conceptually a generalisation of the maximum likelihood ensemble filter by Zupanski (2005).

In an ensemble DA system consisting of $N_e$ ensemble members with an $N$-dimensional state vector, the ensemble mean of the posterior analysis $\bar{x}^a \in \mathbb{R}^N$ is given by $\bar{x}^a = \bar{x}^f + \mathbf{A}w_{\min}$, with $\bar{x}^f \in \mathbb{R}^N$ being the a priori ensemble mean, $\mathbf{A} \in \mathbb{R}^{N \times N_e}$ the matrix of ensemble anomalies with $N$ rows and $N_e$ columns obtained by removing the ensemble mean from the full ensemble matrix, and $w_{\min} \in \mathbb{R}^{N_e}$ the minimum of $w$ obtained from the cost function as

$$w_{\min} := \operatorname*{argmin}_{w} \mathcal{J}(w), \qquad (1a)$$

with

$$\mathcal{J}(w) = \frac{1}{2} (y - \mathcal{H}(\bar{x}^f + \mathbf{A}w))^T \mathbf{R}^{-1} (y - \mathcal{H}(\bar{x}^f + \mathbf{A}w)) + \frac{1}{2}(N-1)w^T w, \qquad (1b)$$

where $y \in \mathbb{R}^{N_0}$ contains $N_0$ observations whose error covariance is specified by $\mathbf{R} \in \mathbb{R}^{N_0 \times N_0}$ and $\mathcal{H}$ is the observation operator. Formulating the problem as in Eq. (1) is equivalent to an ensemble-variational method. Note that extensive reviews of these methods can be found in Chap. 7 of Asch et al. (2016), Bannister (2017) or Sect. 4 in Carrassi et al. (2018). In our implementation, the cost function in Eq. (1) is

minimised using a Gauss–Newton method. The stopping criterion in this study comprises a maximum of 40 iterations for the experiments of state estimation only (Sects. 4.3.2 and 5.1) plus the constraint of $||w^k - w^{k-1}|| < 10^{-3}$ when parameters are also estimated (Sect. 5.2–5.4). The latter additional constraint has been included on the basis of theoretical arguments and numerical experiments. It avoids overfitting at a single analysis step for parameters without time dependency. Without the additional stopping criterion, the parameter estimation leads to excessive corrections that do not appear in the case of state estimation only. The use of early stopping criteria in the context of ensemble-variational methods with domain localisation was also originally suggested by Bocquet (2016) to deal with potential convergence problems. After minimisation, the posterior analysis error covariance is approximated, via the ensemble, by the inverse of the Hessian matrix of $\mathcal{J}$.

The forecast error covariance is approximated with the ensemble anomaly matrix, $\mathbf{A}$, by $\frac{1}{N-1}\mathbf{A}\mathbf{A}^T$. Given that the state vector contains both observed and unobserved model fields, the forecast error covariance matrix contains the (ensemble-based) cross-covariance between these fields that allows for propagating the data content to all model fields, including those that are not directly observed. As we shall clarify later, we use a fully multivariate augmented DA where the state vector contains all model fields and model parameters that will be estimated (see for instance Ruiz et al., 2013). The analysis of the unobserved fields and parameters depends directly on the cross-correlations between these fields and the observations.

## 3 Model setup

In this study, we use a dynamics-only sea-ice model in an idealised setup. In this section, we provide details about the model setup used in our numerical experiments. This includes the model equations, its parameters and the choice of external forcing fields. The setup described in this section is used as modelled “truth” in the following DA experiments.

### 3.1 The dynamics-only sea-ice model

The dynamics-only sea-ice model uses the MEB rheology proposed by Dansereau et al. (2016). The model itself and its numerical implementation, based on an Eulerian, discontinuous Galerkin approach, are presented in Dansereau et al. (2017). In the MEB model, sea ice is treated as a viscous–elasto-brittle material. This rheology allows for representing the ice cover as a brittle solid where it is relatively undamaged (i.e. unfractured) and highly concentrated relative to ice-free waters and as an elastic–viscous fluid where it is intensively fractured and low in concentration. By associating sea-ice with an elastic solid rather than a highly viscous fluid and by incorporating a variable (the level of damage,
both multivariate state estimation in a “perfect-model” setup and parameter estimation in a biased model are reported in Sect. 5. In Sect. 6 we discuss some of the critical aspects of ensemble DA in this context together with how to address them in future works. We finally summarise our findings in Sect. 7.

## 2 Data assimilation: the iterative ensemble Kalman filter

The iterative ensemble Kalman filter (IEnKF; Sakov et al., 2012b; Bocquet and Sakov, 2012, 2014) is a variant of the Kalman filter (KF). Like the KF, it is constructed based on the Bayes theorem under the assumption that the prior (forecast), evidence (observation) and posterior (analysis) follow a Gaussian distribution. It has successfully been applied to joint state and model parameter estimation problems (Bocquet and Sakov, 2013; Haussaire and Bocquet, 2016; Bocquet et al., 2021). In ensemble DA methods, the linear model assumption of the KF is relaxed and the prior distribution is estimated from a finite ensemble of model forecasts. A relevant feature of the IEnKF is that it solves for the analysis via a nonlinear optimisation aimed at maximising the a posteriori probability distribution. This key aspect makes it worth investigating its performance in the context of predicting Arctic sea ice, which is characterised by strong and complex nonlinear relations as well as weak nonlinearity in the observations. The IEnKF is the filter version of the more general iterative ensemble Kalman smoother (Bocquet and Sakov, 2014), and it is conceptually a generalisation of the maximum likelihood ensemble filter by Zupanski (2005).

In an ensemble DA system consisting of $N_e$ ensemble members with an $N$-dimensional state vector, the ensemble mean of the posterior analysis $\bar{x}^a \in \mathbb{R}^N$ is given by $\bar{x}^a = \bar{x}^f + \mathbf{A}w_{\min}$, with $\bar{x}^f \in \mathbb{R}^N$ being the a priori ensemble mean, $\mathbf{A} \in \mathbb{R}^{N \times N_e}$ the matrix of ensemble anomalies with $N$ rows and $N_e$ columns obtained by removing the ensemble mean from the full ensemble matrix, and $w_{\min} \in \mathbb{R}^{N_e}$ the minimum of $w$ obtained from the cost function as

$$w_{\min} := \operatorname{argmin}_{w} \mathcal{J}(w), \qquad (1a)$$

with

$$\mathcal{J}(\mathbf{w}) = \frac{1}{2} (\mathbf{y} - \mathcal{H}(\bar{\mathbf{x}}^f + \mathbf{A}\mathbf{w}))^\mathrm{T} \mathbf{R}^{-1} (\mathbf{y} - \mathcal{H}(\bar{\mathbf{x}}^f + \mathbf{A}\mathbf{w})) \\ + \frac{1}{2}(N-1)\mathbf{w}^\mathrm{T}\mathbf{w}, \qquad (1b)$$

where $\mathbf{y} \in \mathbb{R}^{N_0}$ contains $N_0$ observations whose error covariance is specified by $\mathbf{R} \in \mathbb{R}^{N_0 \times N_0}$ and $\mathcal{H}$ is the observation operator. Formulating the problem as in Eq. (1) is equivalent to an ensemble-variational method. Note that extensive reviews of these methods can be found in Chap. 7 of Asch et al. (2016), Bannister (2017) or Sect. 4 in Carrassi et al. (2018). In our implementation, the cost function in Eq. (1) is

minimised using a Gauss–Newton method. The stopping criterion in this study comprises a maximum of 40 iterations for the experiments of state estimation only (Sects. 4.3.2 and 5.1) plus the constraint of $||\mathbf{w}^k - \mathbf{w}^{k-1}|| < 10^{-3}$ when parameters are also estimated (Sect. 5.2–5.4). The latter additional constraint has been included on the basis of theoretical arguments and numerical experiments. It avoids overfitting at a single analysis step for parameters without time dependency. Without the additional stopping criterion, the parameter estimation leads to excessive corrections that do not appear in the case of state estimation only. The use of early stopping criteria in the context of ensemble-variational methods with domain localisation was also originally suggested by Bocquet (2016) to deal with potential convergence problems. After minimisation, the posterior analysis error covariance is approximated, via the ensemble, by the inverse of the Hessian matrix of $\mathcal{J}$.

The forecast error covariance is approximated with the ensemble anomaly matrix, $\mathbf{A}$, by $\frac{1}{N-1}\mathbf{A}\mathbf{A}^\mathrm{T}$. Given that the state vector contains both observed and unobserved model fields, the forecast error covariance matrix contains the (ensemble-based) cross-covariance between these fields that allows for propagating the data content to all model fields, including those that are not directly observed. As we shall clarify later, we use a fully multivariate augmented DA where the state vector contains all model fields and model parameters that will be estimated (see for instance Ruiz et al., 2013). The analysis of the unobserved fields and parameters depends directly on the cross-correlations between these fields and the observations.

## 3 Model setup

In this study, we use a dynamics-only sea-ice model in an idealised setup. In this section, we provide details about the model setup used in our numerical experiments. This includes the model equations, its parameters and the choice of external forcing fields. The setup described in this section is used as modelled “truth” in the following DA experiments.

### 3.1 The dynamics-only sea-ice model

The dynamics-only sea-ice model uses the MEB rheology proposed by Dansereau et al. (2016). The model itself and its numerical implementation, based on an Eulerian, discontinuous Galerkin approach, are presented in Dansereau et al. (2017). In the MEB model, sea ice is treated as a viscous–elasto-brittle material. This rheology allows for representing the ice cover as a brittle solid where it is relatively undamaged (i.e. unfractured) and highly concentrated relative to ice-free waters and as an elastic–viscous fluid where it is intensively fractured and low in concentration. By associating sea-ice with an elastic solid rather than a highly viscous fluid and by incorporating a variable (the level of damage,
d) to represent its degree of fracturing at the sub-grid scale, this model differs significantly from the widely employed (elastic–)viscous–plastic rheologies (Hunke et al., 2010).

The dynamics-only MEB sea-ice model describes the evolution of nine model fields: sea-ice concentration (SIC) A; sea-ice thickness (SIT) h; sea-ice velocity (SIV) u = (u_x, u_y) = (u, v); level of damage d; cohesion C; and internal stress σ = (σ_{xx} σ_{xy} σ_{yx} σ_{yy}), where σ_{yx} = σ_{xy} (see Dansereau et al., 2017).

In the MEB rheology, the evolution of SIV, the level of damage and stress describes the kinematics of the sea ice. To obtain a physically plausible solution, a stress–velocity– damage constraint that respects the MEB constitutive equa- tion (which relates SIV to σ ) must be satisfied. Other model fields also enter the momentum, constitutive and damage evolution and mass conservation equations, but these are just advected by SIV. Here, we highlight the momentum and stress equations of the sea-ice model. These equations will be relevant in parameter estimation experiments. The momen- tum equation is given as

$$
\rho h \frac{\mathrm{Du}}{Dt} = \nabla \cdot (h\sigma) + \rho_a C_a A |\mathbf{u}_a| \mathbf{u}_a, \quad (2)
$$

where ρ and ρ_a are the sea ice and air density respectively, u_a is the wind field, D/Dt is the material derivative, and C_a is the air drag coefficient. The stress equation is

$$
\lambda^0 d^{\alpha-1} \left[ \frac{\partial \sigma}{\partial t} + (\mathbf{u} \cdot \nabla)\sigma + \beta_a (\nabla \mathbf{u}, \sigma) \right] + [1-d]\sigma = \eta^0 d'^{\alpha} \exp[-c^*(1-A)] \mathbf{K} : \mathcal{D}(\mathbf{u}), \quad (3)
$$

where λ^0 is the undamaged relaxation time; β_a is a function that accounts for the effects of rotation and deformation; η^0 is the undamaged apparent viscosity; d'^α = (1 − η_min/η_0)d^α + η_min/η_0, with η_min being the minimum apparent viscosity and α being a damage-related parameter; K is a stiffness tensor; and D() is the symmetric part of the velocity gradient.

The model equations are discretised on an unstructured tri-
angular grid using a finite-element, discontinuous Galerkin
method where the sea-ice velocities are defined on triangular
vertices (degree of polynomial approximations of 1) and all
other model fields are defined on the face of the triangular el-
ement (degree of polynomial approximations of 0 or constant
by element; see Dansereau et al., 2017). In particular, be-
cause the constitutive, momentum and level-of-damage evo-
lution is coupled, SIV, the level of damage and sea-ice stress
are solved using a semi-implicit, iterative method. Besides
these dynamical properties, the sea-ice model used in this
study does not include thermodynamics processes. As such,
it is a short-timescale proxy for the dynamic–thermodynamic
sea-ice model currently under construction in the SASIP
project, neXtSIM_DG. However, this future model builds on
the MEB rheology and discontinuous Galerkin-based numer-
ical scheme. Therefore, the present model is a reliable surro-

gate of the dynamical components and of the numerical core
of neXtSIM<sub>DG</sub>.

The evolution of the sea-ice model is controlled by multi-
ple model parameters which define the physical properties of
the sea ice and thereby intrinsically affect the representation
of sea ice in the model. These parameters (see Table A1) fol-
low the choice in Dansereau et al. (2017) with exceptions for
the spatial and temporal resolutions and minimum cohesion.
In our experiments, the model runs at the spatial resolution of
around 15 km and a time step of 30 s. This ensures numerical
stability while sufficiently resolving the propagation of dam-
age, the fastest process represented in the model. The model
is solved on a squared model domain of [-L, L] x [-L, L]
with the x and y axes aligned along the perpendicular sides of
the square with L = 100 km. This gives us 512 elements and
285 nodes. To enhance internal variability (and thus maintain
ensemble spread in the IEnKF), the minimum sea-ice cohe-
sion, which sets the resistance of the sea-ice cover to fractur-
ing, is lowered to 5000 from 8000 Pa, which has been used
so far in neXtSIM (Rampal et al., 2016).

In this study, the DA’s ability to estimate two model pa-
rameters will be investigated: the air drag coefficient, C_a,
in Eq. (2) and damage parameter, α, in Eq. (3). C_a enters
the model in the momentum equation (see Eq. 2), modulat-
ing the influence of the wind fields on SIV. This parameter
corresponds to the effect of external forcing on the sea-ice
model. α controls the swiftness of the transition between the
elasto-brittle regime at a low level of damage and the vis-
cous regime at a high level of damage. The mechanical be-
haviour of an elasto-brittle solid is less predictable in nature
than that of a viscous fluid (Weiss and Dansereau, 2017).
This means that erroneous α changes the internal property of
sea ice and influences the sea-ice predictability. The choice
of the parameters to be estimated is largely based on pre-
vious studies. Massonnet et al. (2014) showed that a cou-
pled ocean–sea-ice model is sensitive to the drag coefficient,
which is estimated with the ensemble Kalman filter. Using
neXtSIM, a rheology-based model similar to the MEB model
used here, Rabatel et al. (2018) showed the importance of
the wind stress, and thus the air drag coefficient, in deter-
mining the spread of the sea-ice trajectories. Furthermore,
Cheng et al. (2020) demonstrated that the shape and orienta-
tion of the area covered by the sea-ice trajectories depend on
the sea-ice air drag coefficient. Besides the drag coefficient,
the damage parameter is the only other minimal parameter
in a solid-like MEB rheology-based model. Drag coefficient
and damage parameters are the two sole tunable parameters
present in the MEB model used in this study. This is because,
in a dynamics-only MEB model, the air drag coefficient and
the damage parameter are the only two parameters that do not
affect the maximum speed of the fastest-propagating elastic
waves which influence the model stability with given time
steps.

We have performed a numerical sensitivity analysis which
shows that the observations are sensitive to both parame-
d) to represent its degree of fracturing at the sub-grid scale, this model differs significantly from the widely employed (elastic–)viscous–plastic rheologies (Hunke et al., 2010).

The dynamics-only MEB sea-ice model describes the evolution of nine model fields: sea-ice concentration (SIC) A; sea-ice thickness (SIT) h; sea-ice velocity (SIV) u = (u_x, u_y) = (u, v); level of damage d; cohesion C; and internal stress σ = (σ_{xx} σ_{xy} σ_{yx} σ_{yy}), where σ_{yx} = σ_{xy} (see Dansereau et al., 2017).

In the MEB rheology, the evolution of SIV, the level of damage and stress describes the kinematics of the sea ice. To obtain a physically plausible solution, a stress–velocity– damage constraint that respects the MEB constitutive equa- tion (which relates SIV to σ ) must be satisfied. Other model fields also enter the momentum, constitutive and damage evolution and mass conservation equations, but these are just advected by SIV. Here, we highlight the momentum and stress equations of the sea-ice model. These equations will be relevant in parameter estimation experiments. The momen- tum equation is given as

$$
\rho h \frac{\mathrm{Du}}{Dt} = \nabla \cdot (h\sigma) + \rho_a C_a A |\mathbf{u}_a| \mathbf{u}_a, \quad (2)
$$

where ρ and ρ_a are the sea ice and air density respectively, u_a is the wind field, D/Dt is the material derivative, and C_a is the air drag coefficient. The stress equation is

$$
\lambda^0 d^{\alpha-1} \left[ \frac{\partial \sigma}{\partial t} + (\mathbf{u} \cdot \nabla)\sigma + \beta_a (\nabla \mathbf{u}, \sigma) \right] + [1-d]\sigma = \eta^0 d^{\alpha'} \exp[-c^*(1-A)] \mathbf{K} : \mathcal{D}(\mathbf{u}), \quad (3)
$$

where λ^0 is the undamaged relaxation time; β_a is a function that accounts for the effects of rotation and deformation; η^0 is the undamaged apparent viscosity; d'^α = (1 − η_min/η_0)d^α + η_min/η_0, with η_min being the minimum apparent viscosity and α being a damage-related parameter; K is a stiffness tensor; and D() is the symmetric part of the velocity gradient.

The model equations are discretised on an unstructured tri-
angular grid using a finite-element, discontinuous Galerkin
method where the sea-ice velocities are defined on triangular
vertices (degree of polynomial approximations of 1) and all
other model fields are defined on the face of the triangular el-
ement (degree of polynomial approximations of 0 or constant
by element; see Dansereau et al., 2017). In particular, be-
cause the constitutive, momentum and level-of-damage evo-
lution is coupled, SIV, the level of damage and sea-ice stress
are solved using a semi-implicit, iterative method. Besides
these dynamical properties, the sea-ice model used in this
study does not include thermodynamics processes. As such,
it is a short-timescale proxy for the dynamic–thermodynamic
sea-ice model currently under construction in the SASIP
project, neXtSIM_DG. However, this future model builds on
the MEB rheology and discontinuous Galerkin-based numer-
ical scheme. Therefore, the present model is a reliable surro-

gate of the dynamical components and of the numerical core
of neXtSIM<sub>DG</sub>.

The evolution of the sea-ice model is controlled by multi-
ple model parameters which define the physical properties of
the sea ice and thereby intrinsically affect the representation
of sea ice in the model. These parameters (see Table A1) fol-
low the choice in Dansereau et al. (2017) with exceptions for
the spatial and temporal resolutions and minimum cohesion.
In our experiments, the model runs at the spatial resolution of
around 15 km and a time step of 30 s. This ensures numerical
stability while sufficiently resolving the propagation of dam-
age, the fastest process represented in the model. The model
is solved on a squared model domain of [-L, L] x [-L, L]
with the x and y axes aligned along the perpendicular sides of
the square with L = 100 km. This gives us 512 elements and
285 nodes. To enhance internal variability (and thus maintain
ensemble spread in the IEnKF), the minimum sea-ice cohe-
sion, which sets the resistance of the sea-ice cover to fractur-
ing, is lowered to 5000 from 8000 Pa, which has been used
so far in neXtSIM (Rampal et al., 2016).

In this study, the DA’s ability to estimate two model pa-
rameters will be investigated: the air drag coefficient, C_a,
in Eq. (2) and damage parameter, α, in Eq. (3). C_a enters
the model in the momentum equation (see Eq. 2), modulat-
ing the influence of the wind fields on SIV. This parameter
corresponds to the effect of external forcing on the sea-ice
model. α controls the swiftness of the transition between the
elasto-brittle regime at a low level of damage and the vis-
cous regime at a high level of damage. The mechanical be-
haviour of an elasto-brittle solid is less predictable in nature
than that of a viscous fluid (Weiss and Dansereau, 2017).
This means that erroneous α changes the internal property of
sea ice and influences the sea-ice predictability. The choice
of the parameters to be estimated is largely based on pre-
vious studies. Massonnet et al. (2014) showed that a cou-
pled ocean–sea-ice model is sensitive to the drag coefficient,
which is estimated with the ensemble Kalman filter. Using
neXtSIM, a rheology-based model similar to the MEB model
used here, Rabatel et al. (2018) showed the importance of
the wind stress, and thus the air drag coefficient, in deter-
mining the spread of the sea-ice trajectories. Furthermore,
Cheng et al. (2020) demonstrated that the shape and orienta-
tion of the area covered by the sea-ice trajectories depend on
the sea-ice air drag coefficient. Besides the drag coefficient,
the damage parameter is the only other minimal parameter
in a solid-like MEB rheology-based model. Drag coefficient
and damage parameters are the two sole tunable parameters
present in the MEB model used in this study. This is because,
in a dynamics-only MEB model, the air drag coefficient and
the damage parameter are the only two parameters that do not
affect the maximum speed of the fastest-propagating elastic
waves which influence the model stability with given time
steps.

We have performed a numerical sensitivity analysis which
shows that the observations are sensitive to both parame-
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Sensitivity of SIV, SIC and SIT to perturbed $C_a$ or $\alpha$ parameter. The shaded area is the standard deviation, a representation of uncertainty in the ensemble an indication of the strength of the sensitivity.</figcaption></figure>

ters. Results are shown in Fig. 1, where the same initial and boundary conditions are used while each ensemble member has different $C_a$ or $\alpha$ values. The parameter values are sampled from $N(2.5 \times 10^{-3}, 5 \times 10^{-4})$, and $N(6, 1.5)$ for $C_a$ and $\alpha$ respectively. The model parameters are time-independent, and the ensemble means are $2.5 \times 10^{-3}$ and 6.5 respectively. Figure 1 shows that both $C_a$ and $\alpha$ have a strong impact on SIC and SIT, while a smaller effect is observed on SIV.

The initial and boundary conditions of the model are shown in Table 1. The simulations start with a domain covered by undamaged sea ice at rest with a random cohesion field. The value of the cohesion field for each element is sampled from a uniform distribution between 5000 and 10000 Pa. As shown in Fig. 2e, the initial SIT ($h$) is a “blob” defined by a cosine function given in Table 1, which represents the naturally inhomogeneous distribution of the SIT in space. We use no-slip boundary conditions at $x = -L = L$ and Neumann boundary conditions at $y = -L = L$, where 1 m thick undamaged sea ice is transported into the domain based on SIV. To avoid an influx of sea ice with a uniform cohesion field from the domain boundaries, the cohesion field of the inflowing ice is randomly sampled from a uniform distribution between 5000 and 10000 Pa.

## 3.2 The external wind field

In the model setup, following Dansereau et al. (2017), the ocean is assumed to be at rest and the sea-ice variability relies on the wind forcing. We design an external forcing in the

form of a prescribed storm-like wind drag that mimics the wind field encountered in operational sea-ice DA. Consistent with reality, the wind is therefore a major source of variability in the present simulations (Guemas et al., 2016; Rabatel et al., 2018).

The storm-like wind field is inspired by the test cases for linear advection schemes in Lauritzen et al. (2012) generated by analytical formulae as outlined in Appendix B. We deem this forcing an adequate representation of the complex horizontal atmospheric flow. The wind field is updated at each model computational time step such that the sea ice is always driven by the up-to-date wind field similarly to the case where the sea-ice model is coupled to an atmospheric model.

The cyclone, shown in Fig. 2a, covers around a quarter of the domain in space. The wind field is superposed on a background flow ($2 \text{ m s}^{-1}$) from the bottom to the top of the domain. The centre of the storm moves from $y=0$ to $y=L$ with a speed of $0.25 \text{ m s}^{-1}$, which is slower than the background flow. The parameters of the storm can be found in Table 2. The formulae for the wind field also allow for a fine control of the generation and dissipation of the cyclone. Therefore, as shown in Fig. 2b, the storms (appearing as peaks in the time series) are not persistent but have time-evolving features. Following the dissipation of a cyclone, a “peaceful” period with only the background wind is used. Hence, during our 90 d experiment period, there are 12 storm occurrences. Moreover, like the duration, strength, initial centre position and travel speed of the storm are speci-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sensitivity of SIV, SIC and SIT to perturbed $C_a$ or $\alpha$ parameter. The shaded area is the standard deviation, a representation of uncertainty in the ensemble an indication of the strength of the sensitivity.</figcaption></figure>

ters. Results are shown in Fig. 1, where the same initial and boundary conditions are used while each ensemble member has different $C_a$ or $\alpha$ values. The parameter values are sampled from $N(2.5 \times 10^{-3}, 5 \times 10^{-4})$, and $N(6, 1.5)$ for $C_a$ and $\alpha$ respectively. The model parameters are time-independent, and the ensemble means are $2.5 \times 10^{-3}$ and 6.5 respectively. Figure 1 shows that both $C_a$ and $\alpha$ have a strong impact on SIC and SIT, while a smaller effect is observed on SIV.

The initial and boundary conditions of the model are shown in Table 1. The simulations start with a domain covered by undamaged sea ice at rest with a random cohesion field. The value of the cohesion field for each element is sampled from a uniform distribution between 5000 and 10000 Pa. As shown in Fig. 2e, the initial SIT ($h$) is a “blob” defined by a cosine function given in Table 1, which represents the naturally inhomogeneous distribution of the SIT in space. We use no-slip boundary conditions at $x = -L = L$ and Neumann boundary conditions at $y = -L = L$, where 1 m thick undamaged sea ice is transported into the domain based on SIV. To avoid an influx of sea ice with a uniform cohesion field from the domain boundaries, the cohesion field of the inflowing ice is randomly sampled from a uniform distribution between 5000 and 10000 Pa.

## 3.2 The external wind field

In the model setup, following Dansereau et al. (2017), the ocean is assumed to be at rest and the sea-ice variability relies on the wind forcing. We design an external forcing in the

form of a prescribed storm-like wind drag that mimics the wind field encountered in operational sea-ice DA. Consistent with reality, the wind is therefore a major source of variability in the present simulations (Guemas et al., 2016; Rabatel et al., 2018).

The storm-like wind field is inspired by the test cases for linear advection schemes in Lauritzen et al. (2012) generated by analytical formulae as outlined in Appendix B. We deem this forcing an adequate representation of the complex horizontal atmospheric flow. The wind field is updated at each model computational time step such that the sea ice is always driven by the up-to-date wind field similarly to the case where the sea-ice model is coupled to an atmospheric model.

The cyclone, shown in Fig. 2a, covers around a quarter of the domain in space. The wind field is superposed on a background flow ($2 \text{ m s}^{-1}$) from the bottom to the top of the domain. The centre of the storm moves from $y=0$ to $y=L$ with a speed of $0.25 \text{ m s}^{-1}$, which is slower than the background flow. The parameters of the storm can be found in Table 2. The formulae for the wind field also allow for a fine control of the generation and dissipation of the cyclone. Therefore, as shown in Fig. 2b, the storms (appearing as peaks in the time series) are not persistent but have time-evolving features. Following the dissipation of a cyclone, a “peaceful” period with only the background wind is used. Hence, during our 90 d experiment period, there are 12 storm occurrences. Moreover, like the duration, strength, initial centre position and travel speed of the storm are speci-
<table><caption>**Table 1.** The initial and boundary condition of the experimental setup in the dynamics-only sea-ice model. Here $r = \sqrt{x^2 + y^2}$, with $(x, y)$ being the coordinate of the grid points and $r_0 = 50$ km. The boundary condition at $y = -L = L$ is that of the fields transported into the model domain.</caption><thead><tr><th rowspan="2">Initial condition</th><th colspan="2">Boundary condition</th><th colspan="2">Forcing</th></tr><tr><th>x = -L = L</th><th>y = -L = L</th><th>Ocean</th><th>Atmosphere</th></tr></thead><tbody><tr><td>A = 1</td><td></td><td>A = 1</td><td></td><td></td></tr><tr><td>d = 0</td><td></td><td>d = 0</td><td></td><td></td></tr><tr><td><b>u</b> = 0 m s<sup>-1</sup></td><td rowspan="4"><b>u</b> = 0 m s<sup>-1</sup></td><td>σ(t) · <b>n</b> = 0 Pa</td><td rowspan="4">at rest</td><td>random quasi-periodic storm-like</td></tr><tr><td>σ = 0 Pa</td><td>σ = 0 Pa</td><td></td></tr><tr><td>h = max(1, 1 + cos(π/2 r/r<sub>0</sub>)) m</td><td>h = 1 m</td><td>wind field occurrence</td></tr><tr><td>C ~ <b>U</b>(5000, 10 000) Pa</td><td>C ~ <b>U</b>(5000, 10 000) Pa</td><td></td></tr></tbody></table>

<table><caption>**Table 2.** The parameters of the storm that served as the truth of the wind field.</caption><thead><tr><th rowspan="2">Period (days)</th><th colspan="3">Strength (m s<sup>-1</sup>)</th><th colspan="2">Initial centre (km)</th><th rowspan="2">Travel speed (m s<sup>-1</sup>)</th></tr><tr><th>Rotational</th><th>Divergent</th><th>Background</th><th>x</th><th>y</th></tr></thead><tbody><tr><td>U(2, 5)</td><td>22</td><td>0.1</td><td>2</td><td>U(-25, 25)</td><td>U(-50, 20)</td><td>0.25</td></tr></tbody></table>

fied by the analytical formula with storm duration and initial position sampled from a uniform distribution as shown in Table 2. With the given random initial position of the storm, the variability in the wind field is confined in a limited region of the domain as shown in Fig. 2c and d. Figure 2e–f show that the variability in the wind field leads to variability in the modelled sea ice, mainly in the upper-right region of the domain as the sea ice is mostly damaged in these regions.

## 4 Data assimilation setup

The operational sea-ice DA corrects only the model state without dealing with model errors explicitly. It is, however, impossible to obtain a perfect model. Model errors originate from diverse sources, such as numerical discretisation, the lack of sufficient resolution as well as errors in the model parameters. We focus here on the parametric errors, as these are particularly problematic and dominate the initial conditions' errors in long-term forecasts.

As introduced in Sect. 2, the IEnKF can infer unobserved model fields thanks to its ensemble-based cross-correlations with the observed fields. This feature can also be exploited for parameter estimations. In the parameter estimation, we use the augmented approach in which the model parameters are included as part of the state vector. For example, in the case of estimating the air drag coefficient $C_a$, the state vector is $\mathbf{x} = (h, u, v, \dots, C_a)^T$. Thus their inference is ultimately related to the capability of the IEnKF to properly describe the correlation within this new, augmented virtual state composed by the physical variables and the parameters.

The simultaneous state and parameter estimation is more challenging than state estimation alone, since changes in the

parameters have a direct impact and can substantially modify the model's dynamical properties. In the worst-case scenario, inappropriate parameter values can push the model through a bifurcation (a tipping point in the case of non-autonomous systems), leading to qualitatively different model behaviour than that of the data. Moreover, different parameters can have aliased uncertainties in the observed model fields, implying that not all of the parameters can be uniquely identified.

These challenges motivate the investigation of state and parameter estimation in the idealised dynamics-only sea-ice model. In this study, a series of twin experiments are conducted whereby a model run is taken to represent the truth. Synthetic observations are generated from the truth following a specified observation error distribution.

To assess DA's ability, four different scenarios are explored: (1) a perfect model where the parameters are equal to their “true” values; (2) a model with a biased air drag coefficient, $C_a$; (3) a model with a biased damage parameter, $\alpha$; and (4) a model with biased $C_a$ and $\alpha$. A summary of the setup of each experiment is given in Table 3.

### 4.1 Ensemble generation

The IEnKF belongs to the category of the ensemble-based DA methods. As such it relies on an ensemble of model trajectories to approximate the forecast uncertainty. The ensemble spread represents the error in the estimate of the model state and parameters. We will explore four different scenarios, whereby, as we shall clarify later, we will employ a different strategy each time to generate the ensemble. In particular, in the cases with parametric error, each member of the ensemble will be given a different set of model parameters.
<table><caption>**Table 1.** The initial and boundary condition of the experimental setup in the dynamics-only sea-ice model. Here $r = \sqrt{x^2 + y^2}$, with $(x, y)$ being the coordinate of the grid points and $r_0 = 50$ km. The boundary condition at $y = -L = L$ is that of the fields transported into the model domain.</caption><thead><tr><th rowspan="2">Initial condition</th><th colspan="2">Boundary condition</th><th colspan="2">Forcing</th></tr><tr><th>x = -L = L</th><th>y = -L = L</th><th>Ocean</th><th>Atmosphere</th></tr></thead><tbody><tr><td>A = 1</td><td></td><td>A = 1</td><td></td><td></td></tr><tr><td>d = 0</td><td></td><td>d = 0</td><td></td><td></td></tr><tr><td><b>u</b> = 0 m s<sup>-1</sup></td><td rowspan="4"><b>u</b> = 0 m s<sup>-1</sup></td><td>σ(t) · <b>n</b> = 0 Pa</td><td rowspan="4">at rest</td><td>random quasi-periodic storm-like</td></tr><tr><td>σ = 0 Pa</td><td>σ = 0 Pa</td><td></td></tr><tr><td>h = max(1, 1 + cos(π/2 r/r<sub>0</sub>)) m</td><td>h = 1 m</td><td>wind field occurrence</td></tr><tr><td>C ~ <b>U</b>(5000, 10 000) Pa</td><td>C ~ <b>U</b>(5000, 10 000) Pa</td><td></td></tr></tbody></table>

<table><caption>**Table 2.** The parameters of the storm that served as the truth of the wind field.</caption><thead><tr><th rowspan="2">Period (days)</th><th colspan="3">Strength (m s<sup>-1</sup>)</th><th colspan="2">Initial centre (km)</th><th rowspan="2">Travel speed (m s<sup>-1</sup>)</th></tr><tr><th>Rotational</th><th>Divergent</th><th>Background</th><th>x</th><th>y</th></tr></thead><tbody><tr><td>U(2, 5)</td><td>22</td><td>0.1</td><td>2</td><td>U(-25, 25)</td><td>U(-50, 20)</td><td>0.25</td></tr></tbody></table>

fied by the analytical formula with storm duration and initial position sampled from a uniform distribution as shown in Table 2. With the given random initial position of the storm, the variability in the wind field is confined in a limited region of the domain as shown in Fig. 2c and d. Figure 2e–f show that the variability in the wind field leads to variability in the modelled sea ice, mainly in the upper-right region of the domain as the sea ice is mostly damaged in these regions.

## 4 Data assimilation setup

The operational sea-ice DA corrects only the model state without dealing with model errors explicitly. It is, however, impossible to obtain a perfect model. Model errors originate from diverse sources, such as numerical discretisation, the lack of sufficient resolution as well as errors in the model parameters. We focus here on the parametric errors, as these are particularly problematic and dominate the initial conditions' errors in long-term forecasts.

As introduced in Sect. 2, the IEnKF can infer unobserved model fields thanks to its ensemble-based cross-correlations with the observed fields. This feature can also be exploited for parameter estimations. In the parameter estimation, we use the augmented approach in which the model parameters are included as part of the state vector. For example, in the case of estimating the air drag coefficient $C_a$, the state vector is $\mathbf{x} = (h, u, v, \dots, C_a)^T$. Thus their inference is ultimately related to the capability of the IEnKF to properly describe the correlation within this new, augmented virtual state composed by the physical variables and the parameters.

The simultaneous state and parameter estimation is more challenging than state estimation alone, since changes in the

parameters have a direct impact and can substantially modify the model's dynamical properties. In the worst-case scenario, inappropriate parameter values can push the model through a bifurcation (a tipping point in the case of non-autonomous systems), leading to qualitatively different model behaviour than that of the data. Moreover, different parameters can have aliased uncertainties in the observed model fields, implying that not all of the parameters can be uniquely identified.

These challenges motivate the investigation of state and parameter estimation in the idealised dynamics-only sea-ice model. In this study, a series of twin experiments are conducted whereby a model run is taken to represent the truth. Synthetic observations are generated from the truth following a specified observation error distribution.

To assess DA's ability, four different scenarios are explored: (1) a perfect model where the parameters are equal to their “true” values; (2) a model with a biased air drag coefficient, $C_a$; (3) a model with a biased damage parameter, $\alpha$; and (4) a model with biased $C_a$ and $\alpha$. A summary of the setup of each experiment is given in Table 3.

### 4.1 Ensemble generation

The IEnKF belongs to the category of the ensemble-based DA methods. As such it relies on an ensemble of model trajectories to approximate the forecast uncertainty. The ensemble spread represents the error in the estimate of the model state and parameters. We will explore four different scenarios, whereby, as we shall clarify later, we will employ a different strategy each time to generate the ensemble. In particular, in the cases with parametric error, each member of the ensemble will be given a different set of model parameters.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** An illustration of the wind field and the model response in SIT. (a) A snapshot of the wind field at an arbitrary time; (b) time series of the spatially averaged wind velocity; (c) standard deviation of the *u* component and (d) the *v* component of the wind in time over 90 d; (e) the initial condition of the SIT; (f) SIT after 42 d of simulation.</figcaption></figure>

<table><caption>**Table 3.** The experiment setup used in four different scenarios. These experiments use the same truth, where **u**<sub>a</sub> is the wind forcing and **u**(*t* = 0) is the initial condition of SIV.</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Model error</th><th rowspan="2">Inference vector</th><th rowspan="2">Perturbations</th><th colspan="2">Background</th><th colspan="3">Localisation (km)</th></tr><tr><th>C<sub>a</sub></th><th>α</th><th>SIC</th><th>SIT</th><th>SIV</th></tr></thead><tbody><tr><td>1</td><td>none</td><td>state</td><td>u, u<sub>a</sub></td><td>1.5 × 10<sup>-3</sup></td><td>4</td><td>5</td><td>200 and 30</td><td>30</td></tr><tr><td>2</td><td>C<sub>a</sub></td><td>state, C<sub>a</sub></td><td>u, u<sub>a</sub>, C<sub>a</sub></td><td>2.5 × 10<sup>-3</sup></td><td>4</td><td>5</td><td>30</td><td>30</td></tr><tr><td>3</td><td>α</td><td>state, α</td><td>u, u<sub>a</sub>, α</td><td>1.5 × 10<sup>-3</sup></td><td>6.5</td><td>5</td><td>30</td><td>30</td></tr><tr><td>4</td><td>C<sub>a</sub> and α</td><td>state, C<sub>a</sub>, α</td><td>u, u<sub>a</sub>, C<sub>a</sub>, α</td><td>2.5 × 10<sup>-3</sup></td><td>6.5</td><td>5</td><td>30</td><td>30</td></tr></tbody></table>

Nevertheless, in all of the four scenarios, we will perturb the wind field (i.e. an external forcing), the initial condition of the sea-ice velocity and the cohesion flux at model boundaries. The external atmospheric wind forcing is a major source of forecast uncertainty in sea-ice models. For example, Rabatel et al. (2018) and Cheng et al. (2020) studied the

sensitivity of neXtSIM to the wind forcing and sea-ice cohesion. In this study, we generate synthetic perturbed wind fields around the “true wind” defined in Sect. 3.2 and use them to form an ensemble. The duration, strength, initial position and travel speed of each cyclone are perturbed with noise sampled from Gaussian distributions (see Table 4 for
<figure><img src="image_5.png" /><figcaption>**Figure 2.** An illustration of the wind field and the model response in SIT. (a) A snapshot of the wind field at an arbitrary time; (b) time series of the spatially averaged wind velocity; (c) standard deviation of the *u* component and (d) the *v* component of the wind in time over 90 d; (e) the initial condition of the SIT; (f) SIT after 42 d of simulation.</figcaption></figure>

<table><caption>**Table 3.** The experiment setup used in four different scenarios. These experiments use the same truth, where **u**<sub>a</sub> is the wind forcing and **u**(*t* = 0) is the initial condition of SIV.</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Model error</th><th rowspan="2">Inference vector</th><th rowspan="2">Perturbations</th><th colspan="2">Background</th><th colspan="3">Localisation (km)</th></tr><tr><th>C<sub>a</sub></th><th>α</th><th>SIC</th><th>SIT</th><th>SIV</th></tr></thead><tbody><tr><td>1</td><td>none</td><td>state</td><td>u, u<sub>a</sub></td><td>1.5 × 10<sup>-3</sup></td><td>4</td><td>5</td><td>200 and 30</td><td>30</td></tr><tr><td>2</td><td>C<sub>a</sub></td><td>state, C<sub>a</sub></td><td>u, u<sub>a</sub>, C<sub>a</sub></td><td>2.5 × 10<sup>-3</sup></td><td>4</td><td>5</td><td>30</td><td>30</td></tr><tr><td>3</td><td>α</td><td>state, α</td><td>u, u<sub>a</sub>, α</td><td>1.5 × 10<sup>-3</sup></td><td>6.5</td><td>5</td><td>30</td><td>30</td></tr><tr><td>4</td><td>C<sub>a</sub> and α</td><td>state, C<sub>a</sub>, α</td><td>u, u<sub>a</sub>, C<sub>a</sub>, α</td><td>2.5 × 10<sup>-3</sup></td><td>6.5</td><td>5</td><td>30</td><td>30</td></tr></tbody></table>

Nevertheless, in all of the four scenarios, we will perturb the wind field (i.e. an external forcing), the initial condition of the sea-ice velocity and the cohesion flux at model boundaries. The external atmospheric wind forcing is a major source of forecast uncertainty in sea-ice models. For example, Rabatel et al. (2018) and Cheng et al. (2020) studied the

sensitivity of neXtSIM to the wind forcing and sea-ice cohesion. In this study, we generate synthetic perturbed wind fields around the “true wind” defined in Sect. 3.2 and use them to form an ensemble. The duration, strength, initial position and travel speed of each cyclone are perturbed with noise sampled from Gaussian distributions (see Table 4 for
details). In addition to these perturbations, we also introduce a random walk for the centre of the cyclone in the zonal direction as part of the ensemble perturbation. This is generated as red noise at each time step according to

$$ \Delta x_{i+1} = e^{-\Delta t/\tau} \Delta x_i + \epsilon, \quad (4) $$

where $\Delta x_i$ is the distance travelled zonally at the ith time step, $\tau = 60$ is a time decorrelation factor and $\epsilon$ is noise sampled from $N(0, 4.44 \times 10^{-5})$ m. This red noise is not applied to the true wind. To avoid the cyclone travelling outside of our model domain, we resample the noise if the centre of the cyclone goes outside of the region between -40 and 50 km for the x coordinate. The same check and resampling are applied to the perturbations for the initial centre of the cyclone.

Albeit marginally with respect to variability due to the external forcing, the internal model variability caused by nonlinearities is another source of forecast error. This is accounted for in our ensemble by perturbing the SIV initial condition, in addition to perturbing the external atmospheric wind field. SIV perturbations are sampled from the Gaussian distribution, $N(0, 0.05)$. Furthermore, as discussed in Sect. 3.1, the boundary condition of the random cohesion influx, which differs for each ensemble member, adds another source of uncertainty in our experiments, but as the sea ice drifts slowly and does not travel long distances across the domain, the impact of the cohesion perturbation is limited.

In the experiments with parametric error in either $C_a$ or $\alpha$ or both (experiments 2–4 in Table 3), the prior covariances of the model parameters need to be specified. As shown in Table 4, the parameter values used by the ensemble members are sampled from zero-mean Gaussian distributions with standard deviations set to be around 33% of the true parameter value. Given that both $C_a$ and $\alpha$ are bounded from below due to physical constraints, the sampled values of $C_a$ and $\alpha$ are ensured to be greater than $10^{-5}$ and 2 by rejecting outliers.

## 4.2 Synthetic observations

Synthetic observations in our twin experiments are generated by sampling from the truth with the aim of mimicking how observations of the Arctic sea-ice are collected operationally, i.e. their spatio-temporal density.

Satellite observations of SIC, SIT and sea-ice drift are used in sea-ice forecasting and reanalysis systems (Sakov et al., 2012a). Most operational systems assimilate gridded products. However, some recent studies show the possibility of directly assimilating along-track instead of gridded data for SIT (Fiedler et al., 2022a), which could allow for a more frequent SIT assimilation than the 7 d averaged gridded products (Ricker et al., 2017) because gridded data require time for collection and processing. Hence, a star-shaped spatial distribution mimicking the along-track SIT observations is used here. For simplicity, we assume that the same star-shaped spatial distribution is available daily due to the time

<figure><img src="image_3.png" /><figcaption>Figure 3. Observation distribution of SIT and SIV with the background triangles being the model grid. The observation distribution of SIC is not shown because it covers the entire domain.</figcaption></figure>

required for data collection from the polar satellite. Note that this treatment neglects the temporal variability in the satellite tracks. Following the protocol of gridded data, we generate synthetic observations of SIC and SIV quasi-uniformly across the model domain.

Nowadays SIC satellite observations reach a resolution as high as 10 km. Considering that our model has a spatial resolution of around 15 km, we synthetically observe every grid point for SIC in our experiments. SIV data are sparser, with a spatial resolution of around 50 km. The position of the along-track SIT data is parametrised by four lines on the domain satisfying the condition $|y_i + cx_i| \le r$, where $(x_i, y_i)$ is the spatial coordinate of the grid points, and the pair $(c, r) \in \{(2.2, 16 km), (-2.2, 16 km), (0.5, 8 km), (-0.5, 8 km)\}$. A graphic distribution of SIC and SIV can be found in Fig. 3.

The synthetic observations are generated by sampling from the truth and adding an observational error drawn from Gaussian distributions. The observation error variances of these Gaussians follow those of realistic DA systems. For SIC, we adopt the observation error standard deviation from Sakov et al. (2012a) and Cheng et al. (2023):

$$ \sigma_A[m] = \sqrt{0.01 + (0.5 - |0.5 - A|)^2}. \quad (5) $$

The observation error in the along-track SIT follows the formula for the measurement uncertainty for CryoSat-2 used by Fiedler et al. (2022a):

$$ \sigma_h = \begin{cases} 8 & h < 0.7 \text{ m} \\ 2 \left(1 - \frac{7e}{0.3 - h^e}\right) \frac{h}{100} & 0.7 \le h < 3 \text{ m} \\ \left(5(h-3) + \left(1 - \frac{7e}{0.3 - 3^e}\right)\right) \frac{h}{100} & h > 3 \text{ m.} \end{cases} \quad (6) $$

This equation parametrises the uncertainty in the observations based on the measured value of the SIT. The use of
details). In addition to these perturbations, we also introduce a random walk for the centre of the cyclone in the zonal direction as part of the ensemble perturbation. This is generated as red noise at each time step according to

$$ \Delta x_{i+1} = e^{-\Delta t/\tau} \Delta x_i + \epsilon, \quad (4) $$

where $\Delta x_i$ is the distance travelled zonally at the ith time step, $\tau = 60$ is a time decorrelation factor and $\epsilon$ is noise sampled from $N(0, 4.44 \times 10^{-5})$ m. This red noise is not applied to the true wind. To avoid the cyclone travelling outside of our model domain, we resample the noise if the centre of the cyclone goes outside of the region between -40 and 50 km for the x coordinate. The same check and resampling are applied to the perturbations for the initial centre of the cyclone.

Albeit marginally with respect to variability due to the external forcing, the internal model variability caused by nonlinearities is another source of forecast error. This is accounted for in our ensemble by perturbing the SIV initial condition, in addition to perturbing the external atmospheric wind field. SIV perturbations are sampled from the Gaussian distribution, $N(0, 0.05)$. Furthermore, as discussed in Sect. 3.1, the boundary condition of the random cohesion influx, which differs for each ensemble member, adds another source of uncertainty in our experiments, but as the sea ice drifts slowly and does not travel long distances across the domain, the impact of the cohesion perturbation is limited.

In the experiments with parametric error in either $C_a$ or $\alpha$ or both (experiments 2–4 in Table 3), the prior covariances of the model parameters need to be specified. As shown in Table 4, the parameter values used by the ensemble members are sampled from zero-mean Gaussian distributions with standard deviations set to be around 33% of the true parameter value. Given that both $C_a$ and $\alpha$ are bounded from below due to physical constraints, the sampled values of $C_a$ and $\alpha$ are ensured to be greater than $10^{-5}$ and 2 by rejecting outliers.

## 4.2 Synthetic observations

Synthetic observations in our twin experiments are generated by sampling from the truth with the aim of mimicking how observations of the Arctic sea-ice are collected operationally, i.e. their spatio-temporal density.

Satellite observations of SIC, SIT and sea-ice drift are used in sea-ice forecasting and reanalysis systems (Sakov et al., 2012a). Most operational systems assimilate gridded products. However, some recent studies show the possibility of directly assimilating along-track instead of gridded data for SIT (Fiedler et al., 2022a), which could allow for a more frequent SIT assimilation than the 7 d averaged gridded products (Ricker et al., 2017) because gridded data require time for collection and processing. Hence, a star-shaped spatial distribution mimicking the along-track SIT observations is used here. For simplicity, we assume that the same star-shaped spatial distribution is available daily due to the time

<figure><img src="image_6.png" /><figcaption>Figure 3. Observation distribution of SIT and SIV with the background triangles being the model grid. The observation distribution of SIC is not shown because it covers the entire domain.</figcaption></figure>

required for data collection from the polar satellite. Note that this treatment neglects the temporal variability in the satellite tracks. Following the protocol of gridded data, we generate synthetic observations of SIC and SIV quasi-uniformly across the model domain.

Nowadays SIC satellite observations reach a resolution as high as 10 km. Considering that our model has a spatial resolution of around 15 km, we synthetically observe every grid point for SIC in our experiments. SIV data are sparser, with a spatial resolution of around 50 km. The position of the along-track SIT data is parametrised by four lines on the domain satisfying the condition $|y_i + cx_i| \le r$, where $(x_i, y_i)$ is the spatial coordinate of the grid points, and the pair $(c, r) \in \{(2.2, 16 km), (-2.2, 16 km), (0.5, 8 km), (-0.5, 8 km)\}$. A graphic distribution of SIC and SIV can be found in Fig. 3.

The synthetic observations are generated by sampling from the truth and adding an observational error drawn from Gaussian distributions. The observation error variances of these Gaussians follow those of realistic DA systems. For SIC, we adopt the observation error standard deviation from Sakov et al. (2012a) and Cheng et al. (2023):

$$ \sigma_A[m] = \sqrt{0.01 + (0.5 - |0.5 - A|)^2}. \quad (5) $$

The observation error in the along-track SIT follows the formula for the measurement uncertainty for CryoSat-2 used by Fiedler et al. (2022a):

$$ \sigma_h = \begin{cases} 8 & h < 0.7 \text{ m} \\ 2 \left(1 - \frac{7e}{0.3 - h^e}\right) \frac{h}{100} & 0.7 \le h < 3 \text{ m} \\ \left(5(h-3) + \left(1 - \frac{7e}{0.3 - 3^e}\right)\right) \frac{h}{100} & h > 3 \text{ m.} \end{cases} \quad (6) $$

This equation parametrises the uncertainty in the observations based on the measured value of the SIT. The use of
<table><caption>**Table 4.** Gaussian distributions of wind fields, SIV and parameter perturbations. The mean of the wind fields is given in Table 2.</caption><thead><tr><th rowspan="3">Period (days)</th><th colspan="4">Wind field</th><th rowspan="3">Travel speed (m s<sup>-1</sup>)</th><th rowspan="3">SIV (m s<sup>-1</sup>)</th><th rowspan="3">C<sub>a</sub></th><th rowspan="3">&alpha;</th></tr><tr><th colspan="2">Strength (m s<sup>-1</sup>)</th><th colspan="2">Initial centre (km)</th></tr><tr><th>Rotational</th><th>Background</th><th>x</th><th>y</th></tr></thead><tbody><tr><td>N(0, 0.5)</td><td>N(0, 0.5)</td><td>N(0, &radic;0.05)</td><td>N(0, 1.5)</td><td>N(0, 1.5)</td><td>N(0, 0.0045)</td><td>N(0, 0.05)</td><td>N(0, 5 &times; 10<sup>-4</sup>)</td><td>N(0, 1.5)</td></tr></tbody></table>

<table><caption>**Table 5.** The time- and space-averaged ensemble standard deviation of each free-run ensemble over 90 d.</caption><thead><tr><th>Scenario</th><th>u (m s<sup>-1</sup>)</th><th>v (m s<sup>-1</sup>)</th><th>A</th><th>h (m)</th><th>d</th><th>C (Pa)</th><th>&sigma;<sub>xx</sub> (Pa)</th><th>&sigma;<sub>xy</sub> (Pa)</th><th>&sigma;<sub>yy</sub> (Pa)</th></tr></thead><tbody><tr><td>1</td><td>7.52 &times; 10<sup>-3</sup></td><td>0.011</td><td>0.072</td><td>0.037</td><td>0.215</td><td>153.128</td><td>2584.221</td><td>1475.218</td><td>1701.683</td></tr><tr><td>2</td><td>0.030</td><td>0.031</td><td>0.153</td><td>0.089</td><td>0.240</td><td>289.626</td><td>3089.592</td><td>1508.822</td><td>1925.960</td></tr><tr><td>3</td><td>7.42 &times; 10<sup>-3</sup></td><td>0.011</td><td>0.074</td><td>0.038</td><td>0.185</td><td>166.199</td><td>2657.060</td><td>1502.896</td><td>1754.829</td></tr><tr><td>4</td><td>0.033</td><td>0.034</td><td>0.149</td><td>0.093</td><td>0.210</td><td>295.344</td><td>3139.044</td><td>1517.597</td><td>1951.430</td></tr></tbody></table>

$\sigma_h = 8$ m for $h < 0.7$ m effectively eliminates SIT observations of thin ice from the DA, reflecting the fact that SIT of thin ice from satellite altimetry is notoriously untrustworthy. The error in Eq. (6) does not account for the representation errors. To account for it, we have therefore added an extra factor of 2.

For the sake of simplicity, synthetic SIV observations are used as a substitute for the sea-ice drift data. There is no explicit equation for the observational error variance of SIV, nor are there prototypical examples due to the lack of investigation of SIV uncertainties. Hence, for each observation grid point, we use either 80% of the standard deviation of a single 90d model trajectory (in our case, the truth) of the observation variance or $8 \times 10^{-4}$ m s$^{-1}$ (~0.21 km per 3d), whichever is greater. This leads to a maximum observational error of around 12km per 3d (13km per 3d) for the $u$ ($v$) component and ensures that SIV observations will contribute to the DA correction. However, the value for the standard deviation is smaller than the 14km per 3d value used in Sakov et al. (2012a): one of the few cases of operational reanalyses using SIV we are aware of. Since Sakov et al. (2012a) showed only a faint impact of SIV assimilation, a reduction in the observation error is expected to take better advantage of the SIV data.

## 4.3 Inflation and localisation

### 4.3.1 Inflation

The application of ensemble DA in geosciences is plagued by sampling errors that arise because of the impossibility of using sufficiently large ensembles. The huge size of realistic numerical models of geofluids and the computational constraints imply that the number of an affordable ensemble size is much smaller than the state vector's dimension, $N_e \ll N$. Inflation and localisation are the two main approaches to alleviate sampling errors (and to some extent, model errors, as discussed in Scheffler et al., 2022, and Grudzien et al., 2018).

Following Cheng et al. (2023), an ensemble size of $N_e = 40$ is used for the IEnKF. The ensemble size is thus orders of magnitude smaller than the size of the state vector, which is $O(10^3)$.

In our four experimental scenarios, we use the adaptive inflation method proposed by Bocquet and Sakov (2012). The IEnKF-N is an extension of the IEnKF that includes an adaptive inflation method designed to counteract the sampling error by keeping a safe ensemble spread. To achieve this, the IEnKF-N introduces an uninformative hyperprior in the EnKF: the formulation has been proved to be equivalent to the multiplicative inflation. Although the IEnKF-N has a slightly larger computational cost, it spares us from the very costly offline tuning of the inflation factor (a procedure that should ideally be repeated for each experimental setup).

### 4.3.2 Localisation

We implement domain localisation as described in Bocquet (2016, his Table 2). In the domain localisation, each model grid point (local domain) assimilates observations within a circle centred on the grid point. To gradually reduce the impact of observations away from the central grid point and to ensure spatially smooth DA corrections, the observation error is tapered by the Gaspari–Cohn (GC) function with a cutoff value of $10^{-5}$. In this way, each local domain has a different cost function based on different observations. Moreover, the localisation radii are dependent on the observation type and the spatial correlation scale of the model physics. The latter is naturally time- and space-dependent. In practice, for the sake of computational efficiency and for the difficulties inherent in its adequate assessment, the localisation radius is often a fixed value (although dependence on the physical variable is usually accommodated).

We anticipate that the localisation radius for global parameters (drag coefficient, damage parameter, etc.) is set to infinity. This choice reflects the fact that the parameter does not have a spatial decorrelation scale and is indeed global.
<table><caption>**Table 4.** Gaussian distributions of wind fields, SIV and parameter perturbations. The mean of the wind fields is given in Table 2.</caption><thead><tr><th rowspan="3">Period (days)</th><th colspan="4">Wind field</th><th rowspan="3">Travel speed (m s<sup>-1</sup>)</th><th rowspan="3">SIV (m s<sup>-1</sup>)</th><th rowspan="3">C<sub>a</sub></th><th rowspan="3">&alpha;</th></tr><tr><th colspan="2">Strength (m s<sup>-1</sup>)</th><th colspan="2">Initial centre (km)</th></tr><tr><th>Rotational</th><th>Background</th><th>x</th><th>y</th></tr></thead><tbody><tr><td>N(0, 0.5)</td><td>N(0, 0.5)</td><td>N(0, &radic;0.05)</td><td>N(0, 1.5)</td><td>N(0, 1.5)</td><td>N(0, 0.0045)</td><td>N(0, 0.05)</td><td>N(0, 5 &times; 10<sup>-4</sup>)</td><td>N(0, 1.5)</td></tr></tbody></table>

<table><caption>**Table 5.** The time- and space-averaged ensemble standard deviation of each free-run ensemble over 90 d.</caption><thead><tr><th>Scenario</th><th>u (m s<sup>-1</sup>)</th><th>v (m s<sup>-1</sup>)</th><th>A</th><th>h (m)</th><th>d</th><th>C (Pa)</th><th>&sigma;<sub>xx</sub> (Pa)</th><th>&sigma;<sub>xy</sub> (Pa)</th><th>&sigma;<sub>yy</sub> (Pa)</th></tr></thead><tbody><tr><td>1</td><td>7.52 &times; 10<sup>-3</sup></td><td>0.011</td><td>0.072</td><td>0.037</td><td>0.215</td><td>153.128</td><td>2584.221</td><td>1475.218</td><td>1701.683</td></tr><tr><td>2</td><td>0.030</td><td>0.031</td><td>0.153</td><td>0.089</td><td>0.240</td><td>289.626</td><td>3089.592</td><td>1508.822</td><td>1925.960</td></tr><tr><td>3</td><td>7.42 &times; 10<sup>-3</sup></td><td>0.011</td><td>0.074</td><td>0.038</td><td>0.185</td><td>166.199</td><td>2657.060</td><td>1502.896</td><td>1754.829</td></tr><tr><td>4</td><td>0.033</td><td>0.034</td><td>0.149</td><td>0.093</td><td>0.210</td><td>295.344</td><td>3139.044</td><td>1517.597</td><td>1951.430</td></tr></tbody></table>

$\sigma_h = 8$ m for $h < 0.7$ m effectively eliminates SIT observations of thin ice from the DA, reflecting the fact that SIT of thin ice from satellite altimetry is notoriously untrustworthy. The error in Eq. (6) does not account for the representation errors. To account for it, we have therefore added an extra factor of 2.

For the sake of simplicity, synthetic SIV observations are used as a substitute for the sea-ice drift data. There is no explicit equation for the observational error variance of SIV, nor are there prototypical examples due to the lack of investigation of SIV uncertainties. Hence, for each observation grid point, we use either 80% of the standard deviation of a single 90d model trajectory (in our case, the truth) of the observation variance or $8 \times 10^{-4}$ m s$^{-1}$ (~0.21 km per 3d), whichever is greater. This leads to a maximum observational error of around 12km per 3d (13km per 3d) for the $u$ ($v$) component and ensures that SIV observations will contribute to the DA correction. However, the value for the standard deviation is smaller than the 14km per 3d value used in Sakov et al. (2012a): one of the few cases of operational reanalyses using SIV we are aware of. Since Sakov et al. (2012a) showed only a faint impact of SIV assimilation, a reduction in the observation error is expected to take better advantage of the SIV data.

## 4.3 Inflation and localisation

### 4.3.1 Inflation

The application of ensemble DA in geosciences is plagued by sampling errors that arise because of the impossibility of using sufficiently large ensembles. The huge size of realistic numerical models of geofluids and the computational constraints imply that the number of an affordable ensemble size is much smaller than the state vector's dimension, $N_e \ll N$. Inflation and localisation are the two main approaches to alleviate sampling errors (and to some extent, model errors, as discussed in Scheffler et al., 2022, and Grudzien et al., 2018).

Following Cheng et al. (2023), an ensemble size of $N_e = 40$ is used for the IEnKF. The ensemble size is thus orders of magnitude smaller than the size of the state vector, which is $O(10^3)$.

In our four experimental scenarios, we use the adaptive inflation method proposed by Bocquet and Sakov (2012). The IEnKF-N is an extension of the IEnKF that includes an adaptive inflation method designed to counteract the sampling error by keeping a safe ensemble spread. To achieve this, the IEnKF-N introduces an uninformative hyperprior in the EnKF: the formulation has been proved to be equivalent to the multiplicative inflation. Although the IEnKF-N has a slightly larger computational cost, it spares us from the very costly offline tuning of the inflation factor (a procedure that should ideally be repeated for each experimental setup).

### 4.3.2 Localisation

We implement domain localisation as described in Bocquet (2016, his Table 2). In the domain localisation, each model grid point (local domain) assimilates observations within a circle centred on the grid point. To gradually reduce the impact of observations away from the central grid point and to ensure spatially smooth DA corrections, the observation error is tapered by the Gaspari–Cohn (GC) function with a cutoff value of $10^{-5}$. In this way, each local domain has a different cost function based on different observations. Moreover, the localisation radii are dependent on the observation type and the spatial correlation scale of the model physics. The latter is naturally time- and space-dependent. In practice, for the sake of computational efficiency and for the difficulties inherent in its adequate assessment, the localisation radius is often a fixed value (although dependence on the physical variable is usually accommodated).

We anticipate that the localisation radius for global parameters (drag coefficient, damage parameter, etc.) is set to infinity. This choice reflects the fact that the parameter does not have a spatial decorrelation scale and is indeed global.
This treatment of localisation for global parameters is also adopted in Aksoy et al. (2006), Ruiz et al. (2013) and Massonnet et al. (2014). See also Ruckstuhl and Janjić (2018), Bocquet et al. (2021), and Malartic et al. (2022) for more recent developments on this topic.

As a trade-off between efficiency and accuracy, we opted for optimising the localisation radius when only one model field is observed. Furthermore, we assume that SIV has the same correlation length scale along both components (i.e. no preferred direction of movement), which is consistent with the experimental setup.

For each of the localisation radii explored, we run a 15 d long simulation after a 42 d ensemble free run without DA. Although 15 d may appear too short for adequate tuning and may not include all possible physical regimes, it still covers multiple storm events (storm passage every 3.5 d on average).

The ratio between the root-mean-square error (RMSE) of the analysis and of the free run, as a function of the localisation radius when assimilating SIT, SIC or SIV, is shown in Fig. 4. When assimilating SIT (Fig. 4a), the analysis error in SIT decreases monotonically with the localisation radius until the localisation radius reaches the domain size (200 km). The inset shows that assimilating SIT also leads to a reduction in SIV's RMSE (RMSEª / RMSEfree systematically below 1), while it leads to a deterioration in SIC as soon as the localisation radius is larger than 5 km.

As opposed to when SIT is observed, the localisation radius for SIC observations (Fig. 4b) is very small. After 5 km, the analysis error in SIC increases monotonically with the localisation radius. Observing SIC improves both SIT and, for very long radii, SIV (see inset in Fig. 4b).

Finally, from Fig. 4c we see that the lowest analysis RMSE when assimilating SIV is attained with a localisation radius of around 30 km and that both SIC and SIT will improve in the multivariate update.

Based on these results, the most effective localisation radii for the observed model fields are 200 km for the SIT observations, 5 km for the SIC observations and 30 km for the SIV observations. These different localisation radii arise from the physical spatial correlations and the observation density of these model fields. A summary of the choices for the localisation radii in this study is given in Table 3, and further discussion will be presented in Sect. 5.1.

## 4.4 Treatment of bounded physical variables and model parameters

Of the nine MEB model variables, three are bounded: SIC, the level of damage and SIT. Inferring them via the IEnKF is therefore challenging given the Gaussian assumption on which the IEnKF based. In practice, there is no guarantee that the update analyses of SIC, damage level or SIT fall within their bounds.

Hence, the solution to that problem that we adopt in this study, although sub-optimal, is very pragmatic and straight-

forward. When the analysis of these variables falls outside of their bounds, they are forcibly set to their nearest bounds. We are conscious that this approach leads to local ensemble collapse (whereby some members originally having out-of-range analysis values are all made equal to the nearest boundary value) and can cause biases in the analysis. Nevertheless, the results in Sect. 5 will demonstrate that, when the bounds are not exceeded too often, the approach works well, it does not cause major ensemble collapse and it is successful in removing nonphysical values.

We follow a similar strategy when performing the estimation of the damage parameter, $\alpha$. There the analysed value of $\alpha$ is assumed to be bounded from below by $\alpha = 2$, the value at which the analysed $\alpha$ will be taken if lower than 2.

Our pragmatic approach is however insufficient when estimating the drag coefficient $C_a$. This parameter is bounded to be strictly positive; therefore ensemble collapse could happen whenever the ensemble mean of $C_a$ approaches zero. In that case a potentially large number of members may receive negative analysis values that would all be restored to the same small positive values. Hence, we adopt a resampling approach in which, with each negative analysis of $C_a$, we sample from $N(10^{-5}, 10^{-6})$ until a positive value is obtained.

## 5 Results

We evaluate the performance of the IEnKF for state and parameter estimation in the dynamics-only MEB sea-ice model under the four different scenarios described in Table 3. We use the root-mean-square error (RMSE) over both time- and space as a skill metric.

### 5.1 Scenario 1: inferring the model's physical variables under a perfect-model setup

Here we study in detail the fully multivariate DA using different combinations of observations under a perfect-model scenario. The RMSE is calculated over 30 d long assimilation experiments that follow a 42 d free ensemble run without DA. In the 30 d assimilation, we assimilate observations daily as mentioned in Sect. 4.2, leading to a total of 31 analyses.

As mentioned in Sect. 4.4, three out of the nine model fields are bounded quantities: SIC, SIT and the level of damage. Given that respect of those bounds is not automatically guaranteed by the DA procedure, we apply a post-processing step. Here, we quantify how often nonphysical values (i.e. values out of the bounds) are produced in the analysis: Table 6 shows the number of violations during the 30 d DA period.

When only SIC is assimilated, physical bounds for SIC and for the level of damage are exceeded very occasionally, for below 1.5% of the total analyses and grid points. On the other hand, when SIC is not observed, the chance of
This treatment of localisation for global parameters is also adopted in Aksoy et al. (2006), Ruiz et al. (2013) and Massonnet et al. (2014). See also Ruckstuhl and Janjić (2018), Bocquet et al. (2021), and Malartic et al. (2022) for more recent developments on this topic.

As a trade-off between efficiency and accuracy, we opted for optimising the localisation radius when only one model field is observed. Furthermore, we assume that SIV has the same correlation length scale along both components (i.e. no preferred direction of movement), which is consistent with the experimental setup.

For each of the localisation radii explored, we run a 15 d long simulation after a 42 d ensemble free run without DA. Although 15 d may appear too short for adequate tuning and may not include all possible physical regimes, it still covers multiple storm events (storm passage every 3.5 d on average).

The ratio between the root-mean-square error (RMSE) of the analysis and of the free run, as a function of the localisation radius when assimilating SIT, SIC or SIV, is shown in Fig. 4. When assimilating SIT (Fig. 4a), the analysis error in SIT decreases monotonically with the localisation radius until the localisation radius reaches the domain size (200 km). The inset shows that assimilating SIT also leads to a reduction in SIV's RMSE (RMSEª / RMSEfree systematically below 1), while it leads to a deterioration in SIC as soon as the localisation radius is larger than 5 km.

As opposed to when SIT is observed, the localisation radius for SIC observations (Fig. 4b) is very small. After 5 km, the analysis error in SIC increases monotonically with the localisation radius. Observing SIC improves both SIT and, for very long radii, SIV (see inset in Fig. 4b).

Finally, from Fig. 4c we see that the lowest analysis RMSE when assimilating SIV is attained with a localisation radius of around 30 km and that both SIC and SIT will improve in the multivariate update.

Based on these results, the most effective localisation radii for the observed model fields are 200 km for the SIT observations, 5 km for the SIC observations and 30 km for the SIV observations. These different localisation radii arise from the physical spatial correlations and the observation density of these model fields. A summary of the choices for the localisation radii in this study is given in Table 3, and further discussion will be presented in Sect. 5.1.

## 4.4 Treatment of bounded physical variables and model parameters

Of the nine MEB model variables, three are bounded: SIC, the level of damage and SIT. Inferring them via the IEnKF is therefore challenging given the Gaussian assumption on which the IEnKF based. In practice, there is no guarantee that the update analyses of SIC, damage level or SIT fall within their bounds.

Hence, the solution to that problem that we adopt in this study, although sub-optimal, is very pragmatic and straight-

forward. When the analysis of these variables falls outside of their bounds, they are forcibly set to their nearest bounds. We are conscious that this approach leads to local ensemble collapse (whereby some members originally having out-of-range analysis values are all made equal to the nearest boundary value) and can cause biases in the analysis. Nevertheless, the results in Sect. 5 will demonstrate that, when the bounds are not exceeded too often, the approach works well, it does not cause major ensemble collapse and it is successful in removing nonphysical values.

We follow a similar strategy when performing the estimation of the damage parameter, $\alpha$. There the analysed value of $\alpha$ is assumed to be bounded from below by $\alpha = 2$, the value at which the analysed $\alpha$ will be taken if lower than 2.

Our pragmatic approach is however insufficient when estimating the drag coefficient $C_a$. This parameter is bounded to be strictly positive; therefore ensemble collapse could happen whenever the ensemble mean of $C_a$ approaches zero. In that case a potentially large number of members may receive negative analysis values that would all be restored to the same small positive values. Hence, we adopt a resampling approach in which, with each negative analysis of $C_a$, we sample from $N(10^{-5}, 10^{-6})$ until a positive value is obtained.

## 5 Results

We evaluate the performance of the IEnKF for state and parameter estimation in the dynamics-only MEB sea-ice model under the four different scenarios described in Table 3. We use the root-mean-square error (RMSE) over both time- and space as a skill metric.

### 5.1 Scenario 1: inferring the model's physical variables under a perfect-model setup

Here we study in detail the fully multivariate DA using different combinations of observations under a perfect-model scenario. The RMSE is calculated over 30 d long assimilation experiments that follow a 42 d free ensemble run without DA. In the 30 d assimilation, we assimilate observations daily as mentioned in Sect. 4.2, leading to a total of 31 analyses.

As mentioned in Sect. 4.4, three out of the nine model fields are bounded quantities: SIC, SIT and the level of damage. Given that respect of those bounds is not automatically guaranteed by the DA procedure, we apply a post-processing step. Here, we quantify how often nonphysical values (i.e. values out of the bounds) are produced in the analysis: Table 6 shows the number of violations during the 30 d DA period.

When only SIC is assimilated, physical bounds for SIC and for the level of damage are exceeded very occasionally, for below 1.5% of the total analyses and grid points. On the other hand, when SIC is not observed, the chance of
<figure><img src="image_4.png" /><figcaption>**Figure 4.** RMSE<sup>a</sup> / RMSE<sup>free</sup> as a function of the localisation radius, in experiments with observed SIT (a), SIC (b) or SIV (c). The red star indicates the lowest analysis error. The insets show the error in two other unobserved model fields out of the three observed model variables, and the dashed red line indicates the point where RMSE<sup>a</sup> = RMSE<sup>free</sup>.</figcaption></figure>

<table><caption>**Table 6.** The total percentage of local analyses that violate the physical bounds in the 30 d multivariate DA experiments. “SIC+SIT30+SIV” is an experiment where SIT uses a 30 km localisation radius, whereas “SIC+SIT+SIV” uses a localisation radius of 200 km for SIT.</caption><thead><tr><th rowspan="2">Assimilated observations</th><th colspan="4">Physical bounds' violation</th></tr><tr><th>SIC &lt; 0</th><th>SIC &gt; 1</th><th>d &lt; 0</th><th>d &gt; 1</th></tr></thead><tbody><tr><td>SIC</td><td>0.01%</td><td>0.01%</td><td>0.02%</td><td>1.41%</td></tr><tr><td>SIT</td><td>0.25%</td><td>6.95%</td><td>0.34%</td><td>18.66%</td></tr><tr><td>SIV</td><td>0.23%</td><td>1.36%</td><td>0.19%</td><td>11.80%</td></tr><tr><td>SIC + SIT</td><td>0.01%</td><td>3.45%</td><td>0.09%</td><td>20.93%</td></tr><tr><td>SIC + SIV</td><td>0.01%</td><td>0.85%</td><td>0.17%</td><td>12.09%</td></tr><tr><td>SIT + SIV</td><td>0.22%</td><td>5.17%</td><td>0.24%</td><td>18.87%</td></tr><tr><td>SIC + SIT + SIV</td><td>0.01%</td><td>3.75%</td><td>0.17%</td><td>20.64%</td></tr><tr><td>SIC + SIT30 + SIV</td><td>0.01%</td><td>0.98%</td><td>0.14%</td><td>12.97%</td></tr></tbody></table>

The bound SIT > 0 is never violated because SIT is always larger than 1 m in our experiments.

getting nonphysical SIC analyses increases by 2 orders of magnitude, although it remains below 5.2%. Similarly, the SIC observations lower the chance of there being nonphysical damage analyses. Notably, when SIT is observed, it leads to analyses that more often violate the physical bounds, particularly the upper bound for the level of damage. The damage bound violations are more severe than those in SIC, since most of the time the sea ice is undamaged and thus very close to the potentially violable bounds of the model fields. Without the possibility of observing the damage field, the cross-correlations may amplify the analysis increments. This effect can be mitigated by limiting the localisation radius of SIT to 30 km as shown in the row of SIC + SIT30 + SIV. The violation of the physical bounds is efficiently addressed by the post-processing step that brings them within their physical limits as described in Sect. 4.4, and it does not lead to an increased RMSE afterwards. As shown in Fig. 5, the SIC analysis is improved compared to the free run in all cases and the sea-ice damage is the only model field that can be less accurate than the free run.

Figure 5 shows the ratios of analysis/free-run (a) and forecast/free-run (b) RMSE (averaged in space and time)

for different types of observations (y axis) and variables (x axis). Values smaller than 1.0 indicate that DA brings about, in general, an improvement in the state estimate compared to the free run. Comparing panels (a) and (b), we can evaluate how much of the DA update is effective in reducing the forecast error (recall that the analysis cycle is 1 d long). From the figure we immediately notice that when only one observation field is assimilated, that same field is the area where most of the improvement in the analysis occurs. This is consistent with our results in Sect. 4.3.2. Figure 5 also shows almost no changes/slight improvements in SIV analysis when SIC alone is assimilated and conversely for the SIC analysis when SIT is assimilated. This appears to be in contrast to Fig. 4, which indicates that SIT observations have a negative impact on the SIC analysis and that SIC observations can cause SIV analysis to deteriorate at certain localisation radii (see Fig. 4a and b insets). Nevertheless, the longer experiments of Fig. 5 suggest that the ensemble has acquired dynamical consistency and therefore better reproduces the cross-variable correlations.

In the MEB rheology, the level of damage, cohesion and stress have a close relationship. Without assimilating SIT,
<figure><img src="image_7.png" /><figcaption>**Figure 4.** RMSE<sup>a</sup> / RMSE<sup>free</sup> as a function of the localisation radius, in experiments with observed SIT (a), SIC (b) or SIV (c). The red star indicates the lowest analysis error. The insets show the error in two other unobserved model fields out of the three observed model variables, and the dashed red line indicates the point where RMSE<sup>a</sup> = RMSE<sup>free</sup>.</figcaption></figure>

<table><caption>**Table 6.** The total percentage of local analyses that violate the physical bounds in the 30 d multivariate DA experiments. “SIC+SIT30+SIV” is an experiment where SIT uses a 30 km localisation radius, whereas “SIC+SIT+SIV” uses a localisation radius of 200 km for SIT.</caption><thead><tr><th rowspan="2">Assimilated observations</th><th colspan="4">Physical bounds' violation</th></tr><tr><th>SIC &lt; 0</th><th>SIC &gt; 1</th><th>d &lt; 0</th><th>d &gt; 1</th></tr></thead><tbody><tr><td>SIC</td><td>0.01%</td><td>0.01%</td><td>0.02%</td><td>1.41%</td></tr><tr><td>SIT</td><td>0.25%</td><td>6.95%</td><td>0.34%</td><td>18.66%</td></tr><tr><td>SIV</td><td>0.23%</td><td>1.36%</td><td>0.19%</td><td>11.80%</td></tr><tr><td>SIC + SIT</td><td>0.01%</td><td>3.45%</td><td>0.09%</td><td>20.93%</td></tr><tr><td>SIC + SIV</td><td>0.01%</td><td>0.85%</td><td>0.17%</td><td>12.09%</td></tr><tr><td>SIT + SIV</td><td>0.22%</td><td>5.17%</td><td>0.24%</td><td>18.87%</td></tr><tr><td>SIC + SIT + SIV</td><td>0.01%</td><td>3.75%</td><td>0.17%</td><td>20.64%</td></tr><tr><td>SIC + SIT30 + SIV</td><td>0.01%</td><td>0.98%</td><td>0.14%</td><td>12.97%</td></tr></tbody></table>

The bound SIT > 0 is never violated because SIT is always larger than 1 m in our experiments.

getting nonphysical SIC analyses increases by 2 orders of magnitude, although it remains below 5.2%. Similarly, the SIC observations lower the chance of there being nonphysical damage analyses. Notably, when SIT is observed, it leads to analyses that more often violate the physical bounds, particularly the upper bound for the level of damage. The damage bound violations are more severe than those in SIC, since most of the time the sea ice is undamaged and thus very close to the potentially violable bounds of the model fields. Without the possibility of observing the damage field, the cross-correlations may amplify the analysis increments. This effect can be mitigated by limiting the localisation radius of SIT to 30 km as shown in the row of SIC + SIT30 + SIV. The violation of the physical bounds is efficiently addressed by the post-processing step that brings them within their physical limits as described in Sect. 4.4, and it does not lead to an increased RMSE afterwards. As shown in Fig. 5, the SIC analysis is improved compared to the free run in all cases and the sea-ice damage is the only model field that can be less accurate than the free run.

Figure 5 shows the ratios of analysis/free-run (a) and forecast/free-run (b) RMSE (averaged in space and time)

for different types of observations (y axis) and variables (x axis). Values smaller than 1.0 indicate that DA brings about, in general, an improvement in the state estimate compared to the free run. Comparing panels (a) and (b), we can evaluate how much of the DA update is effective in reducing the forecast error (recall that the analysis cycle is 1 d long). From the figure we immediately notice that when only one observation field is assimilated, that same field is the area where most of the improvement in the analysis occurs. This is consistent with our results in Sect. 4.3.2. Figure 5 also shows almost no changes/slight improvements in SIV analysis when SIC alone is assimilated and conversely for the SIC analysis when SIT is assimilated. This appears to be in contrast to Fig. 4, which indicates that SIT observations have a negative impact on the SIC analysis and that SIC observations can cause SIV analysis to deteriorate at certain localisation radii (see Fig. 4a and b insets). Nevertheless, the longer experiments of Fig. 5 suggest that the ensemble has acquired dynamical consistency and therefore better reproduces the cross-variable correlations.

In the MEB rheology, the level of damage, cohesion and stress have a close relationship. Without assimilating SIT,
<figure><img src="image_5.png" /><figcaption>**Figure 5.** RMSE in the state estimate for scenario 1. (a) RMSEª/ RMSEªfree. (b) RMSEf/ RMSEffree. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** A histogram of the ensemble correlation between model variables when all observations are assimilated over the period of assimilation. The legend indicates the standard deviation of the histogram. A larger spread/standard deviation represents an overall higher correlation between variables.</figcaption></figure>

the IEnKF-N leads to improvements in sea-ice cohesion and stress and some improvements in the level of damage too. However, the assimilation of SIT tends to result in overly damaged sea ice. This is also evident from Table 6, where it can be seen that the assimilation of SIT leads to higher chances of breaking the bounds of the level of damage. This adverse effect on the damage and stress persists when SIT is assimilated together with SIV. Interestingly, however, when SIT is assimilated together with SIC, the boundedness of the level of damage is improved for undamaged sea ice ($d < 0$) but is not for completely damaged sea ice ($d > 1$). However, it is sufficient to improve the overall RMSE of the level of damage (see Fig. 5). One possible reason for this is that, without the thermodynamics, the forecast error mainly comes from the damaged sea ice and the overestimation of undam-

aged sea ice makes little contribution to the RMSE after the post-processing. An ad hoc remedy is to use a shorter localisation radius for SIT (experiment SIC + SIT30 + SIV in Fig. 5), although it also reduces the improvements in SIT and cohesion. Moreover, as a result of better cross-correlations, the use of small localisation mitigates the violations of physical bounds, as shown in Table 6. Another approach to avoiding the negative impact of SIT consists in artificially setting the covariance matrix entries between SIT and SIC and between SIT and damage to zero. This is achieved in practice by assimilating only SIC and SIV for the SIT and damage variables while all observations are assimilated for the rest of the state vector. The experiment “var loc” shows improved results in Fig. 5.

We show the strength of the cross-correlations between different model variables in Fig. 6. The values are taken from the experiment where all observations are assimilated, from all spatial points and all analyses. As expected, the distributions all peak around zero: this is because beyond a certain distance, the correlations are all very small (a fact that is at the basis of the use of localisation), with the larger values concentrated in proximity to the analysis point and populating mainly the tails of the distributions in Fig. 6. The width of the distributions indicates that in many instances the correlations are (in absolute values) as high as 0.5. To provide a quantitative comparison among the distributions' width, we also show the standard deviations. These cross-correlations can be understood from a physical point of view. The cross-correlations between SIV and stress are the strongest as SIV is mainly driven by the external wind field. In addition, SIV, stress and the level of damage are closely coupled processes, so the level of damage and SIV also show strong cross-correlation. The weak cross-correlation between SIV and SIC and between SIV and SIT is a result of the small magnitude of SIV, which transports SIC, SIT and cohesion. Nevertheless, this gives rise to a strong correlation between SIC
<figure><img src="image_8.png" /><figcaption>**Figure 5.** RMSE in the state estimate for scenario 1. (a) RMSEª/RMSEfree. (b) RMSEf/RMSEfree. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** A histogram of the ensemble correlation between model variables when all observations are assimilated over the period of assimilation. The legend indicates the standard deviation of the histogram. A larger spread/standard deviation represents an overall higher correlation between variables.</figcaption></figure>

the IEnKF-N leads to improvements in sea-ice cohesion and stress and some improvements in the level of damage too. However, the assimilation of SIT tends to result in overly damaged sea ice. This is also evident from Table 6, where it can be seen that the assimilation of SIT leads to higher chances of breaking the bounds of the level of damage. This adverse effect on the damage and stress persists when SIT is assimilated together with SIV. Interestingly, however, when SIT is assimilated together with SIC, the boundedness of the level of damage is improved for undamaged sea ice ($d < 0$) but is not for completely damaged sea ice ($d > 1$). However, it is sufficient to improve the overall RMSE of the level of damage (see Fig. 5). One possible reason for this is that, without the thermodynamics, the forecast error mainly comes from the damaged sea ice and the overestimation of undam-

aged sea ice makes little contribution to the RMSE after the post-processing. An ad hoc remedy is to use a shorter localisation radius for SIT (experiment SIC + SIT30 + SIV in Fig. 5), although it also reduces the improvements in SIT and cohesion. Moreover, as a result of better cross-correlations, the use of small localisation mitigates the violations of physical bounds, as shown in Table 6. Another approach to avoiding the negative impact of SIT consists in artificially setting the covariance matrix entries between SIT and SIC and between SIT and damage to zero. This is achieved in practice by assimilating only SIC and SIV for the SIT and damage variables while all observations are assimilated for the rest of the state vector. The experiment “var loc” shows improved results in Fig. 5.

We show the strength of the cross-correlations between different model variables in Fig. 6. The values are taken from the experiment where all observations are assimilated, from all spatial points and all analyses. As expected, the distributions all peak around zero: this is because beyond a certain distance, the correlations are all very small (a fact that is at the basis of the use of localisation), with the larger values concentrated in proximity to the analysis point and populating mainly the tails of the distributions in Fig. 6. The width of the distributions indicates that in many instances the correlations are (in absolute values) as high as 0.5. To provide a quantitative comparison among the distributions' width, we also show the standard deviations. These cross-correlations can be understood from a physical point of view. The cross-correlations between SIV and stress are the strongest as SIV is mainly driven by the external wind field. In addition, SIV, stress and the level of damage are closely coupled processes, so the level of damage and SIV also show strong cross-correlation. The weak cross-correlation between SIV and SIC and between SIV and SIT is a result of the small magnitude of SIV, which transports SIC, SIT and cohesion. Nevertheless, this gives rise to a strong correlation between SIC
and SIT because they are both controlled by the advection processes.

A physical interpretation can also be invoked to explain why the improvements in SIV analysis do not necessarily translate into improvements in the forecast, as shown in Fig. 5. We argue that this is due to the instant injection of error from the wind field after the assimilation of observations. The correction from the assimilation acts as a perturbation to the model, which induces a model adjustment. In contrast to SIV, the increased error in the analysis of the level of damage, when assimilating SIT, is mitigated by the further damage caused by the wind field. On the other hand, the longer timescale of the variability in SIC, SIT and cohesion makes them less sensitive to the model adjustments, and thus better analysis generally yields better sea-ice forecasts (see Fig. 5).

In summary, our experiments show that the IEnKF can improve both the analysis and the forecast of the MEB sea-ice model. We confirm, in line with previous studies, the necessity of controlling the main source of the uncertainty from the external forcing in the sea-ice model. We argue that a correct wind field can counteract the deterioration of the inaccurate SIV forecast. Furthermore, the results also demonstrate the positive impact of a variable-dependent localisation radius in combating both sampling error and nonlinearities. Moreover, even if the IEnKF suffers from sampling error with a large localisation radius, the deterioration of some of the unobserved fields is marginal, with the RMSE being only 6% larger than the free run. With a reduced localisation radius, the IEnKF improves all the unobserved fields.

## 5.2 Scenario 2: inferring the model's physical variables and the drag coefficient $C_a$

We achieved satisfactory performance in the state estimation of the MEB model using the IEnKF-N under the perfect-model assumption. In the experiments described in this section, we assume that $C_a$ is incorrectly specified and attempt to recover the true value using DA, while all other parameters are perfectly known.

The air drag coefficient, $C_a$, controls the degree to which momentum from the wind is transferred to the sea-ice cover. In our model, similarly to most sea-ice models, it is a constant scalar value that modulates the wind drag; see Eq. (2). Equation (2) shows two main sources of uncertainties in the same term: the wind field, $u_a$, and the drag coefficient, $C_a$. The multiplicative role of $C_a$ makes it an amplification factor for the uncertainty that originates from uncertain wind fields. Given that the wind field is the main source of uncertainties in the sea-ice dynamics, the incorrect specification of $C_a$ also affects the predictability of the sea ice. To see this, note the larger free-run ensemble spread in this scenario compared to the perfect-model scenario in Table 5. Note that the difference in the $C_a$ between scenarios 1 and 2 is as small as $10^{-3}$ (see Table 3).

Similarly to the state estimation experiments, the assimilation starts after a 42 d free ensemble run. Figure 7a shows the time series of the analysis of $C_a$, while Fig. 7b shows the analysis spread in $C_a$ as a function of time. We see that after the drastic correction of the initial bias, the ensemble spread stabilises for all experiments. The ensemble spread converges and stabilises to slightly above zero except for the case where only SIT/SIV is assimilated. In those cases, the ensemble spread is larger, which might be a consequence of the sparser observation density. Due to SIV being sensitive to $C_a$, the estimation of $C_a$ is affected by changes in SIV observations. These changes are captured by the adaptive inflation scheme leading to an increased ensemble spread of $C_a$ toward the end of the experiment period. Another remarkable feature is that in all experiments, $C_a$ drops significantly over the first time steps, thereby approaching (but not necessarily converging to) the true value (red line). This is a clear consequence of a strong correlation between the observed fields and the parameter. Besides this, we then observe different converging values and performance depending on the type of observations assimilated. Except for the case where only SIC is assimilated, the analyses underestimate $C_a$ at the end of our experiment time. The smaller increments of $C_a$, when assimilating SIC, are a result of smaller cross-correlation between SIC and $C_a$ in comparison to SIT and $C_a$ as shown in Fig. 8a and c. Similarly to the increments of $C_a$, the SIT experiment also shows increased increments of the observed fields compared to the SIC observation, as in Fig. 8b and d. This suggests that the ensemble spread in SIT is larger than in SIC. Moreover, the cross-correlation of the ensemble anomaly between SIC/SIT and $C_a$ is spatially inhomogeneous (see Fig. 7a and c). This implies that the error is controlled not solely by the global parameter but also by other spatially dependent processes. One possible explanation is the error in the wind fields. As the ensemble error is primarily driven by the error in the wind field scaled by $C_a$ acting as wind forcing, the cross-correlation between the $C_a$ and the observations may be affected by the error in the wind fields. This suggests that while the IEnKF successfully corrects large biases in $C_a$, it may not be able to correct errors of smaller magnitude equally well as the errors can be aliased with the error in the wind field. In the latter case, the estimation of $C_a$ does not necessarily converge to the “true” $C_a$.

The correction of $C_a$ is also effective in the state estimation as shown in Fig. 9. The IEnKF efficiently controls the RMSE for nearly all model fields regardless of the combination of observations compared to the free run that is instead affected by a large bias in $C_a$. Notably, compared to the free run, SIV is improved not only in the analysis, but also in the forecast. This is different from the multivariate update in Sect. 5.1 because the positively biased $C_a$ amplifies the uncertainties from the wind field, which is then reduced by the corrected $C_a$. The time series in Fig. 10c–d show that the corrected $C_a$ significantly reduces the bias in SIV, but the tran-
and SIT because they are both controlled by the advection processes.

A physical interpretation can also be invoked to explain why the improvements in SIV analysis do not necessarily translate into improvements in the forecast, as shown in Fig. 5. We argue that this is due to the instant injection of error from the wind field after the assimilation of observations. The correction from the assimilation acts as a perturbation to the model, which induces a model adjustment. In contrast to SIV, the increased error in the analysis of the level of damage, when assimilating SIT, is mitigated by the further damage caused by the wind field. On the other hand, the longer timescale of the variability in SIC, SIT and cohesion makes them less sensitive to the model adjustments, and thus better analysis generally yields better sea-ice forecasts (see Fig. 5).

In summary, our experiments show that the IEnKF can improve both the analysis and the forecast of the MEB sea-ice model. We confirm, in line with previous studies, the necessity of controlling the main source of the uncertainty from the external forcing in the sea-ice model. We argue that a correct wind field can counteract the deterioration of the inaccurate SIV forecast. Furthermore, the results also demonstrate the positive impact of a variable-dependent localisation radius in combating both sampling error and nonlinearities. Moreover, even if the IEnKF suffers from sampling error with a large localisation radius, the deterioration of some of the unobserved fields is marginal, with the RMSE being only 6% larger than the free run. With a reduced localisation radius, the IEnKF improves all the unobserved fields.

## 5.2 Scenario 2: inferring the model's physical variables and the drag coefficient $C_a$

We achieved satisfactory performance in the state estimation of the MEB model using the IEnKF-N under the perfect-model assumption. In the experiments described in this section, we assume that $C_a$ is incorrectly specified and attempt to recover the true value using DA, while all other parameters are perfectly known.

The air drag coefficient, $C_a$, controls the degree to which momentum from the wind is transferred to the sea-ice cover. In our model, similarly to most sea-ice models, it is a constant scalar value that modulates the wind drag; see Eq. (2). Equation (2) shows two main sources of uncertainties in the same term: the wind field, $u_a$, and the drag coefficient, $C_a$. The multiplicative role of $C_a$ makes it an amplification factor for the uncertainty that originates from uncertain wind fields. Given that the wind field is the main source of uncertainties in the sea-ice dynamics, the incorrect specification of $C_a$ also affects the predictability of the sea ice. To see this, note the larger free-run ensemble spread in this scenario compared to the perfect-model scenario in Table 5. Note that the difference in the $C_a$ between scenarios 1 and 2 is as small as $10^{-3}$ (see Table 3).

Similarly to the state estimation experiments, the assimilation starts after a 42 d free ensemble run. Figure 7a shows the time series of the analysis of $C_a$, while Fig. 7b shows the analysis spread in $C_a$ as a function of time. We see that after the drastic correction of the initial bias, the ensemble spread stabilises for all experiments. The ensemble spread converges and stabilises to slightly above zero except for the case where only SIT/SIV is assimilated. In those cases, the ensemble spread is larger, which might be a consequence of the sparser observation density. Due to SIV being sensitive to $C_a$, the estimation of $C_a$ is affected by changes in SIV observations. These changes are captured by the adaptive inflation scheme leading to an increased ensemble spread of $C_a$ toward the end of the experiment period. Another remarkable feature is that in all experiments, $C_a$ drops significantly over the first time steps, thereby approaching (but not necessarily converging to) the true value (red line). This is a clear consequence of a strong correlation between the observed fields and the parameter. Besides this, we then observe different converging values and performance depending on the type of observations assimilated. Except for the case where only SIC is assimilated, the analyses underestimate $C_a$ at the end of our experiment time. The smaller increments of $C_a$, when assimilating SIC, are a result of smaller cross-correlation between SIC and $C_a$ in comparison to SIT and $C_a$ as shown in Fig. 8a and c. Similarly to the increments of $C_a$, the SIT experiment also shows increased increments of the observed fields compared to the SIC observation, as in Fig. 8b and d. This suggests that the ensemble spread in SIT is larger than in SIC. Moreover, the cross-correlation of the ensemble anomaly between SIC/SIT and $C_a$ is spatially inhomogeneous (see Fig. 7a and c). This implies that the error is controlled not solely by the global parameter but also by other spatially dependent processes. One possible explanation is the error in the wind fields. As the ensemble error is primarily driven by the error in the wind field scaled by $C_a$ acting as wind forcing, the cross-correlation between the $C_a$ and the observations may be affected by the error in the wind fields. This suggests that while the IEnKF successfully corrects large biases in $C_a$, it may not be able to correct errors of smaller magnitude equally well as the errors can be aliased with the error in the wind field. In the latter case, the estimation of $C_a$ does not necessarily converge to the “true” $C_a$.

The correction of $C_a$ is also effective in the state estimation as shown in Fig. 9. The IEnKF efficiently controls the RMSE for nearly all model fields regardless of the combination of observations compared to the free run that is instead affected by a large bias in $C_a$. Notably, compared to the free run, SIV is improved not only in the analysis, but also in the forecast. This is different from the multivariate update in Sect. 5.1 because the positively biased $C_a$ amplifies the uncertainties from the wind field, which is then reduced by the corrected $C_a$. The time series in Fig. 10c–d show that the corrected $C_a$ significantly reduces the bias in SIV, but the tran-
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Time series of (a) the analysis of $C_a$ using different combinations of observations and (b) the analysis ensemble spread of $C_a$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Cross-correlation between the observation and $C_a$: (a) correlation between SIC and $C_a$ at the last analysis; (b) SIC increment. In both cases it refers to the experiment where only SIC is assimilated; panels (c) and (d) are the same as panels (a) and (b) but with SIT in place of SIC in the experiment where only SIT is assimilated.</figcaption></figure>

sient error in the wind forcing still impacts the accuracy of the SIV forecast.

In Fig. 7a, the best skill in estimating $C_a$ is achieved when assimilating SIC alone due to its higher observation density. However, even with fewer observations, SIV achieves good performance due to a close relationship with the wind field and $C_a$ in Eq. (2). Notably, assimilating only SIV greatly re-

duces the RMSE across all model fields except for the sea-ice stress (see the third row in Fig. 9a). The increased RMSE in the sea-ice stress arises because the DA update of the stress sometimes severely violates the constitutive equation. This can potentially lead to unstable model solutions and model crashes. Model crashes can be avoided by using a large number of iterations for the MEB solver during an extended pe-
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Time series of (a) the analysis of $C_a$ using different combinations of observations and (b) the analysis ensemble spread of $C_a$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Cross-correlation between the observation and $C_a$: (a) correlation between SIC and $C_a$ at the last analysis; (b) SIC increment. In both cases it refers to the experiment where only SIC is assimilated; panels (c) and (d) are the same as panels (a) and (b) but with SIT in place of SIC in the experiment where only SIT is assimilated.</figcaption></figure>

sient error in the wind forcing still impacts the accuracy of the SIV forecast.

In Fig. 7a, the best skill in estimating $C_a$ is achieved when assimilating SIC alone due to its higher observation density. However, even with fewer observations, SIV achieves good performance due to a close relationship with the wind field and $C_a$ in Eq. (2). Notably, assimilating only SIV greatly re-

duces the RMSE across all model fields except for the sea-ice stress (see the third row in Fig. 9a). The increased RMSE in the sea-ice stress arises because the DA update of the stress sometimes severely violates the constitutive equation. This can potentially lead to unstable model solutions and model crashes. Model crashes can be avoided by using a large number of iterations for the MEB solver during an extended pe-
<figure><img src="image_9.png" /><figcaption>**Figure 9.** RMSE in the state estimate for scenario 2. (a) RMSEª/ RMSEªfree. (b) RMSEf/ RMSEffree. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Time series of (a) sea-ice stress in the x direction when only SIV is assimilated, (b) sea-ice stress in the x direction when all observed fields are assimilated, and SIV in the (c) u and (d) v component when all observed fields are assimilated. The black points in the time series are the spatially averaged analysis.</figcaption></figure>

riod after the DA step, but the unphysical update still leads to inaccurate forecasts of the sea-ice stress (see the third row in Fig. 9a). The inaccurate stress forecast can be observed in Fig. 10a, which shows a significant underestimation in the stress time series around day 17. Such underestimation does not occur when all observations are assimilated, as shown in Fig. 10b. We found that the erroneous stress underestimation in the time series occurs when the local element of SIC is reduced by the assimilation (not shown) to the point

that it creates a large SIC gradient and increases the elastic behaviour of the sea ice (see the right-hand side of Eq. 3). Since the temporal variation in SIC is small, this issue persists, and a continued decrease in the sea-ice stress along with the model integration is maintained. This incorrect SIC estimate is remedied in the next DA step, where the multivariate DA restores SIC. Hence, assimilating SIC mitigates the unphysical DA update. Arguably, it is unlikely that only SIV is assimilated in real scenarios, yet our results suggest that it is
<figure><img src="image_12.png" /><figcaption>**Figure 9.** RMSE in the state estimate for scenario 2. (a) RMSEª/ RMSEªfree. (b) RMSEf/ RMSEffree. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Time series of (a) sea-ice stress in the x direction when only SIV is assimilated, (b) sea-ice stress in the x direction when all observed fields are assimilated, and SIV in the (c) u and (d) v component when all observed fields are assimilated. The black points in the time series are the spatially averaged analysis.</figcaption></figure>

riod after the DA step, but the unphysical update still leads to inaccurate forecasts of the sea-ice stress (see the third row in Fig. 9a). The inaccurate stress forecast can be observed in Fig. 10a, which shows a significant underestimation in the stress time series around day 17. Such underestimation does not occur when all observations are assimilated, as shown in Fig. 10b. We found that the erroneous stress underestimation in the time series occurs when the local element of SIC is reduced by the assimilation (not shown) to the point

that it creates a large SIC gradient and increases the elastic behaviour of the sea ice (see the right-hand side of Eq. 3). Since the temporal variation in SIC is small, this issue persists, and a continued decrease in the sea-ice stress along with the model integration is maintained. This incorrect SIC estimate is remedied in the next DA step, where the multivariate DA restores SIC. Hence, assimilating SIC mitigates the unphysical DA update. Arguably, it is unlikely that only SIV is assimilated in real scenarios, yet our results suggest that it is
wise to restrict the multivariate update of the sea-ice concentration in this case.

Our experiments show that the IEnKF is able to reduce the bias in $C_a$ based on available sea-ice observations. The improved $C_a$ estimation significantly reduces the error in the model fields. From the momentum equation, Eq. (2), we know that SIV is directly linked to $C_a$. However, our results show that, although assimilating SIV is crucial for estimating $C_a$, SIC and SIT still improve the estimate of $C_a$ and the state. Importantly, when assimilated in conjunction with SIV, SIC and/or SIT mitigates the imbalance of the constitutive equations.

## 5.3 Scenario 3: inferring the model's physical variables and its erroneous damage parameter $\alpha$

While the drag coefficient, $C_a$, is linked to the external forcing, the internal property of the sea ice is largely controlled by the damage parameter $\alpha$. The damage enters the model in the stress equation, Eq. (3). Although $\alpha$ is added to the model in an ad hoc manner, it plays an essential role in the MEB rheology as it sets the rate at which viscosity decreases with increasing level of damage and thereby controls the transition between the elasto-brittle regime at low damage levels and the viscous regime at high damage levels.

Model sensitivity studies by Dansereau (2016) and Weiss and Dansereau (2017) showed that the value of $\alpha$ is critical in determining the macroscopic mechanical behaviour of the model and that a value of $\alpha \ge 4$ leads to complex, sea-ice-compatible, behaviours. In fact, in our experiments, the truth, i.e. the un-biased model, is set to $\alpha = 4$ (see Table 3). In the DA experiments with the estimate of $\alpha$, we mimic an initial biased estimation of the parameter that is in the same range of sea-ice-compatible behaviour: we choose $\alpha = 6.5$ (see Table 3). Our strategy is realistic as the initial “guess” does not cause the model to behave qualitatively differently from the observations when $\alpha \ge 4$. Drastic changes in the dynamical regimes are also challenging for DA, and they influence the error dynamics between the model parameter and the observed fields. Given that the IEnKF updates (in both the state fields and the parameters) are unbounded, we apply the same post-processing to keep the analyses within physically acceptable bounds, as described in Sect. 4.4.

One of the fundamental challenges in estimating $\alpha$ arises from the nonlinear relationship between $\alpha$ and the observed fields. In particular, the parameter is directly related to the stress field, which is not observable. The complex, nonlinear and indirect nature of these relations can lead to inaccurate (finite-)ensemble-based cross-correlations. Another potential challenge comes from the low sensitivity of the model fields to $\alpha$ compared to the wind field and $C_a$ (see Fig. 1). As shown in Table 5, in a 90 d free run, the ensemble spread of observed model fields is only marginally larger than that in the perfect-model scenario.

Despite these obstacles, the IEnKF shows encouraging results in estimating $\alpha$, as shown in Fig. 11a. Assimilating SIC or SIT alone leads to under- and overestimation of $\alpha$ after 30 d. From Fig. 11b, we see that the ensemble spread in SIT is still relatively large at the end of the experiment. This suggests not only that it is only slightly reduced at the analysis steps but also that further adjustment of $\alpha$ beyond the 30th day is still a possibility. The simultaneous assimilation of SIT and SIC leads to a good approximation of the true value of $\alpha = 4$. Similar results are attained whenever SIC is assimilated (see experiment SIC + SIV or SIC + SIT + SIV). This suggests the crucial relationship between $\alpha$ and SIC, which appears to be the key observation for inferring the damage parameter.

Our results also suggest that observations of SIV cannot be used to retrieve $\alpha$ effectively. In all the experiments with SIV observations, the estimated $\alpha$ gradually approaches the truth until day 18, when it then abruptly diverges away from it. An insight into the reasons behind this sudden change is provided in Fig. 12, which shows the spatial distribution (on the model domain) of the cross-correlation between $\alpha$ and either $u$ or $v$, on days 17 and 18. The cross-correlations between SIV and $\alpha$ flip their signs spatially from day 17 to 18. This is related to the uncertainties in the wind field, which dominate the forecast uncertainty in SIV.

Let us illustrate this issue with simple mathematical arguments. We assume that the nonlinear sea-ice model, $M$, can be approximated by its linearisation $M$ within the forecast interval between two successive analyses. The forecast ensemble anomaly (error) of SIV at time step $k$ can be approximated as

$$ \delta \mathbf{u}_k = \mathbf{M}_k^\alpha \delta \alpha_{k-1} + \mathbf{M}_k^{u_a} \delta \mathbf{u}_{a_{k-1}} + \mathbf{M}_k^x \delta \mathbf{x}_{k-1} + \mathcal{O}(2), \quad (7) $$

where $\delta$ represents the deviation from the ensemble mean and $x$ is the state vector except the damage parameter; the superscript means the model sensitivity to the corresponding variables in which $\mathbf{M}_k^\alpha \in \mathbb{R}^{n \times 1}$, $\mathbf{M}_k^{u_a} \in \mathbb{R}^{n \times m}$ and $\mathbf{M}_k^x \in \mathbb{R}^{n \times l}$, with $n$ being the number of SIV values in the controlled vector, $m$ being the number of elements in the wind field vector and $l$ being the number of elements of model state in the controlled vector. $\mathcal{O}(2)$ represents the high-order terms that are greater than second order. In the EnKF (and thus in the IEnKF), the cross-covariance is estimated from the ensemble anomaly of SIV and the perturbations of $\alpha$, which approximate $\mathbb{E}[\delta \mathbf{u}_k^T \delta \alpha_k]$. This cross-covariance between SIV and $\alpha$ is related to the perturbation (error) from the model state, wind and $\alpha$. In our case, whenever $||\mathbf{M}_k^{u_a} \delta \mathbf{u}_a|| \gg \mathbf{M}_k^\alpha \delta \alpha$, the SIV uncertainty is mainly driven by the wind. This can falsely give a strong cross-correlation between SIV and $\alpha$, producing an incorrect $\alpha$ estimate when it is only based on SIV observations. A similar issue was encountered in Simon and Bertino (2012), where the authors found the parameter estimation challenging when the uncertainty in the parameters showed relatively low uncertainties in the observed fields.
wise to restrict the multivariate update of the sea-ice concentration in this case.

Our experiments show that the IEnKF is able to reduce the bias in $C_a$ based on available sea-ice observations. The improved $C_a$ estimation significantly reduces the error in the model fields. From the momentum equation, Eq. (2), we know that SIV is directly linked to $C_a$. However, our results show that, although assimilating SIV is crucial for estimating $C_a$, SIC and SIT still improve the estimate of $C_a$ and the state. Importantly, when assimilated in conjunction with SIV, SIC and/or SIT mitigates the imbalance of the constitutive equations.

## 5.3 Scenario 3: inferring the model's physical variables and its erroneous damage parameter $\alpha$

While the drag coefficient, $C_a$, is linked to the external forcing, the internal property of the sea ice is largely controlled by the damage parameter $\alpha$. The damage enters the model in the stress equation, Eq. (3). Although $\alpha$ is added to the model in an ad hoc manner, it plays an essential role in the MEB rheology as it sets the rate at which viscosity decreases with increasing level of damage and thereby controls the transition between the elasto-brittle regime at low damage levels and the viscous regime at high damage levels.

Model sensitivity studies by Dansereau (2016) and Weiss and Dansereau (2017) showed that the value of $\alpha$ is critical in determining the macroscopic mechanical behaviour of the model and that a value of $\alpha \ge 4$ leads to complex, sea-ice-compatible, behaviours. In fact, in our experiments, the truth, i.e. the un-biased model, is set to $\alpha = 4$ (see Table 3). In the DA experiments with the estimate of $\alpha$, we mimic an initial biased estimation of the parameter that is in the same range of sea-ice-compatible behaviour: we choose $\alpha = 6.5$ (see Table 3). Our strategy is realistic as the initial “guess” does not cause the model to behave qualitatively differently from the observations when $\alpha \ge 4$. Drastic changes in the dynamical regimes are also challenging for DA, and they influence the error dynamics between the model parameter and the observed fields. Given that the IEnKF updates (in both the state fields and the parameters) are unbounded, we apply the same post-processing to keep the analyses within physically acceptable bounds, as described in Sect. 4.4.

One of the fundamental challenges in estimating $\alpha$ arises from the nonlinear relationship between $\alpha$ and the observed fields. In particular, the parameter is directly related to the stress field, which is not observable. The complex, nonlinear and indirect nature of these relations can lead to inaccurate (finite-)ensemble-based cross-correlations. Another potential challenge comes from the low sensitivity of the model fields to $\alpha$ compared to the wind field and $C_a$ (see Fig. 1). As shown in Table 5, in a 90 d free run, the ensemble spread of observed model fields is only marginally larger than that in the perfect-model scenario.

Despite these obstacles, the IEnKF shows encouraging results in estimating $\alpha$, as shown in Fig. 11a. Assimilating SIC or SIT alone leads to under- and overestimation of $\alpha$ after 30 d. From Fig. 11b, we see that the ensemble spread in SIT is still relatively large at the end of the experiment. This suggests not only that it is only slightly reduced at the analysis steps but also that further adjustment of $\alpha$ beyond the 30th day is still a possibility. The simultaneous assimilation of SIT and SIC leads to a good approximation of the true value of $\alpha = 4$. Similar results are attained whenever SIC is assimilated (see experiment SIC + SIV or SIC + SIT + SIV). This suggests the crucial relationship between $\alpha$ and SIC, which appears to be the key observation for inferring the damage parameter.

Our results also suggest that observations of SIV cannot be used to retrieve $\alpha$ effectively. In all the experiments with SIV observations, the estimated $\alpha$ gradually approaches the truth until day 18, when it then abruptly diverges away from it. An insight into the reasons behind this sudden change is provided in Fig. 12, which shows the spatial distribution (on the model domain) of the cross-correlation between $\alpha$ and either $u$ or $v$, on days 17 and 18. The cross-correlations between SIV and $\alpha$ flip their signs spatially from day 17 to 18. This is related to the uncertainties in the wind field, which dominate the forecast uncertainty in SIV.

Let us illustrate this issue with simple mathematical arguments. We assume that the nonlinear sea-ice model, $M$, can be approximated by its linearisation $M$ within the forecast interval between two successive analyses. The forecast ensemble anomaly (error) of SIV at time step $k$ can be approximated as

$$ \delta \mathbf{u}_k = \mathbf{M}_k^\alpha \delta \alpha_{k-1} + \mathbf{M}_k^{u_a} \delta \mathbf{u}_{a_{k-1}} + \mathbf{M}_k^x \delta \mathbf{x}_{k-1} + \mathcal{O}(2), \quad (7) $$

where $\delta$ represents the deviation from the ensemble mean and $x$ is the state vector except the damage parameter; the superscript means the model sensitivity to the corresponding variables in which $\mathbf{M}_k^\alpha \in \mathbb{R}^{n \times 1}$, $\mathbf{M}_k^{u_a} \in \mathbb{R}^{n \times m}$ and $\mathbf{M}_k^x \in \mathbb{R}^{n \times l}$, with $n$ being the number of SIV values in the controlled vector, $m$ being the number of elements in the wind field vector and $l$ being the number of elements of model state in the controlled vector. $\mathcal{O}(2)$ represents the high-order terms that are greater than second order. In the EnKF (and thus in the IEnKF), the cross-covariance is estimated from the ensemble anomaly of SIV and the perturbations of $\alpha$, which approximate $\mathbb{E}[\delta \mathbf{u}_k^T \delta \alpha_k]$. This cross-covariance between SIV and $\alpha$ is related to the perturbation (error) from the model state, wind and $\alpha$. In our case, whenever $||\mathbf{M}_k^{u_a} \delta \mathbf{u}_a|| \gg \mathbf{M}_k^\alpha \delta \alpha$, the SIV uncertainty is mainly driven by the wind. This can falsely give a strong cross-correlation between SIV and $\alpha$, producing an incorrect $\alpha$ estimate when it is only based on SIV observations. A similar issue was encountered in Simon and Bertino (2012), where the authors found the parameter estimation challenging when the uncertainty in the parameters showed relatively low uncertainties in the observed fields.
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Time series of (a) the analysis of $\alpha$ and (b) its ensemble spread.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 12.** Cross-correlation between the (a, b) u and (c, d) v component of SIV and $\alpha$ on day 17 and 18 when only SIV is assimilated.</figcaption></figure>

Based on this argument, we can show that the sign flip of cross-correlations in Fig. 12 is a result of the change in wind field. In Fig. 13, the cross-correlation between the wind field and $\alpha$ shows the periodic northward travel of the storm on the y axis (see Fig. 13a and b) and the rotation of the storm (see Fig. 13c and d). This matches the changes in cross-correlations in Fig. 12. The incorrect estimate of $\alpha$ highlights the challenges when the primary sources of the uncertainties are external instead of being the model parameters. We ob-

serve that this effect also influences the $\alpha$ estimate when SIV is assimilated with the SIT.

Figure 14 shows the improvements in the state estimate. As expected, the analyses of the observed fields are in general more accurate than the free run. Moreover, compared to the perfect-model scenario, the damage field is now improved relative to the free run in all experiments. The stress field is only moderately improved even though it has a direct relationship with $\alpha$ in Eq. (3).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Time series of (a) the analysis of $\alpha$ and (b) its ensemble spread.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Cross-correlation between the (a, b) u and (c, d) v component of SIV and $\alpha$ on day 17 and 18 when only SIV is assimilated.</figcaption></figure>

Based on this argument, we can show that the sign flip of cross-correlations in Fig. 12 is a result of the change in wind field. In Fig. 13, the cross-correlation between the wind field and $\alpha$ shows the periodic northward travel of the storm on the y axis (see Fig. 13a and b) and the rotation of the storm (see Fig. 13c and d). This matches the changes in cross-correlations in Fig. 12. The incorrect estimate of $\alpha$ highlights the challenges when the primary sources of the uncertainties are external instead of being the model parameters. We ob-

serve that this effect also influences the $\alpha$ estimate when SIV is assimilated with the SIT.

Figure 14 shows the improvements in the state estimate. As expected, the analyses of the observed fields are in general more accurate than the free run. Moreover, compared to the perfect-model scenario, the damage field is now improved relative to the free run in all experiments. The stress field is only moderately improved even though it has a direct relationship with $\alpha$ in Eq. (3).
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Same as Fig. 12 but for the components of the wind forcing in place of those of SIV.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 14.** RMSE in the state estimate for scenario 3. (a) RMSE<sup>a</sup> / RMSE<sup>free</sup>. (b) RMSE<sup>f</sup> / RMSE<sup>free</sup>. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Same as Fig. 12 but for the components of the wind forcing in place of those of SIV.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** RMSE in the state estimate for scenario 3. (a) RMSE<sup>a</sup> / RMSE<sup>free</sup>. (b) RMSE<sup>f</sup> / RMSE<sup>free</sup>. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>
Interestingly, without observing SIV, the assimilation of SIC shows a negative impact on the analysis of SIV. This may be a result of the underestimated $\alpha$ in these experiments. With a low value of $\alpha$, the sea ice has a more elastic behaviour. The elastic motion has a short timescale and is sensitive to the perturbation of the sea-ice variables, making the slowly changed SIC observations unreliable. This is consistent with the finding from Weiss and Dansereau (2017) such that, with high $\alpha$ values and the sea ice transitioned from elastic to viscous behaviour, the system becomes more predictable. When SIV is assimilated, the overestimated $\alpha$ leads to only slight error increases in the SIC analysis without hampering the state estimation of other model variables. This suggests that the model state and parameter are adjusted based on the primary source of error, the wind forcing.

Our results demonstrate the possibility of estimating $\alpha$ successfully using as many as possible of the available observations. It is notable that no single type of observation alone can infer $\alpha$ accurately. We also demonstrate that the IEnKF cannot always identify the correct source of error in a complex environment. In our case, the error from the wind field negatively impacts the estimation of $\alpha$. Interestingly, a deteriorated $\alpha$ analysis does not necessarily lead to a deteriorated state estimation. In contrast, the overestimation of $\alpha$ still moderately improves the sea-ice forecast due to improved predictability.

## 5.4 Scenario 4: inferring the model's physical variables and its erroneous $C_a$ and $\alpha$

In the previous sections we demonstrated that, with sufficient observations, the IEnKF can estimate $C_a$ and $\alpha$ when only one of them is erroneous. Considering that both $C_a$ and $\alpha$ are in the closely related equations for SIV and stress (see Eqs. 2 and 3), it is of interest to investigate the possibility of estimating both model parameters simultaneously.

Figure 15b shows the estimated $C_a$ after 30 d of assimilation. All experiments approach the true value of $C_a$ (although starting from an overestimated value). When SIC and SIV are assimilated together, the IEnKF gives the best estimate of $C_a$, while it underestimates $C_a$ when only SIT and SIV are assimilated. The IEnKF also gives an increase in the estimate of $\alpha$ as shown in Fig. 15a. Although this demonstrates the difficulty in estimating both parameters, it is remarkable that it still leads to an improved forecast of model fields as shown in Fig. 16. This is known as a compensating effect (Bocquet, 2012). This apparent contradiction is related to the non-identifiability of the problem whereby the parameter estimation cannot fully reconstruct the true parameters. In practice, this implies that more than one set of parameters can produce the same observed fields. The identifiability issue can also be viewed from a physical standpoint: $C_a$ determines the motion of the sea ice via the wind forcing. With overestimated a priori $C_a$, the wind can cause fast sea-ice motions. $\alpha$ influences the rate at which sea ice transitions

from an elasto-brittle solid behaviour to a viscous fluid behaviour with an increasing level of damage. As discussed in Sect. 5.3, the increased $\alpha$ leads to more viscous sea ice that is more easily subjected to permanent deformations instead of fast transient elastic deformations. In this situation, the IEnKF controls the fast sea-ice motions by removing the elastic sea-ice deformations.

The identifiability issue can also be associated with incorrect cross-correlations between model fields and model parameters from different sources of errors as described in Eq. (7). We observe that assimilation correctly decreases $C_a$, showing a relatively accurate cross-correlation between the observations and $C_a$. This controls the impact of the external wind forcing, as a positively biased $C_a$ can amplify the uncertainty in the wind field, which increases the term $M^{u_a}\delta u_a$ in Eq. (7). As discussed in Sect. 5.3, the outstanding external wind forcing could lead to erroneous estimation of $\alpha$. Hence, after the external uncertainties are controlled, the ensemble should be able to develop reasonable cross-correlations between $\alpha$ and the observations. To test this, we adopt the following strategy in the 30 d assimilation experiments: (1) constraining the external uncertainty by estimating $C_a$ only for 10 d, (2) developing uncertainty from $\alpha$ by model forecasting without assimilation for 10 d and (3) estimating $\alpha$ for the last 10 d. This is labelled "seq" in Fig. 15. Although $C_a$ is still underestimated, we observe a reduction of error in $\alpha$. Also, as other experiments do not show a decreased $\alpha$ after an underestimated $C_a$, this may suggest the importance of step (2), where the uncertainty arising from $\alpha$ developed. We stress here that the 10 d is simply chosen for convenience and is not tuned.

To further consolidate our findings on the performance of the IEnKF in the simultaneous parameter estimations of $C_a$ and $\alpha$, extensive, and computationally demanding, investigations with different model states, forcings and parameters would be ideal. In fact, even though the model is a simplification of a full pan-Arctic sea-ice model, the computational power quickly scales up to beyond our resources.

With the aforementioned computational constraints in mind, we performed three additional experiments using different truths. In addition to the existing experiment (EXP$\alpha^L C_a^L$) where we assimilate all observations to estimate $\alpha$ and $C_a$, we implemented three other experiments where the truth uses (1) $\alpha = 4$ and $C_a = 3.5 \times 10^{-3}$ (EXP$\alpha^L C_a^H$), (2) $\alpha = 7$ and $C_a = 1.5 \times 10^{-3}$ (EXP$\alpha^H C_a^L$), and (3) $\alpha = 7$ and $C_a = 3.5 \times 10^{-3}$ (EXP$\alpha^H C_a^H$). Here, the superscript L and H denote that the truth is lower and higher respectively than the initial guess. We investigate relevant scenarios (e.g. under- or overestimation of the real values) in which the model's qualitative behaviour is of the same sort. The latter specifically concerns the model structural stability, that is to say the fact that the model is not subject to bifurcation of its general behaviour.

When the truth of model parameters is higher than the initial guess, the air drag coefficient is $2\sigma$ above the initial guess
Interestingly, without observing SIV, the assimilation of SIC shows a negative impact on the analysis of SIV. This may be a result of the underestimated $\alpha$ in these experiments. With a low value of $\alpha$, the sea ice has a more elastic behaviour. The elastic motion has a short timescale and is sensitive to the perturbation of the sea-ice variables, making the slowly changed SIC observations unreliable. This is consistent with the finding from Weiss and Dansereau (2017) such that, with high $\alpha$ values and the sea ice transitioned from elastic to viscous behaviour, the system becomes more predictable. When SIV is assimilated, the overestimated $\alpha$ leads to only slight error increases in the SIC analysis without hampering the state estimation of other model variables. This suggests that the model state and parameter are adjusted based on the primary source of error, the wind forcing.

Our results demonstrate the possibility of estimating $\alpha$ successfully using as many as possible of the available observations. It is notable that no single type of observation alone can infer $\alpha$ accurately. We also demonstrate that the IEnKF cannot always identify the correct source of error in a complex environment. In our case, the error from the wind field negatively impacts the estimation of $\alpha$. Interestingly, a deteriorated $\alpha$ analysis does not necessarily lead to a deteriorated state estimation. In contrast, the overestimation of $\alpha$ still moderately improves the sea-ice forecast due to improved predictability.

## 5.4 Scenario 4: inferring the model's physical variables and its erroneous $C_a$ and $\alpha$

In the previous sections we demonstrated that, with sufficient observations, the IEnKF can estimate $C_a$ and $\alpha$ when only one of them is erroneous. Considering that both $C_a$ and $\alpha$ are in the closely related equations for SIV and stress (see Eqs. 2 and 3), it is of interest to investigate the possibility of estimating both model parameters simultaneously.

Figure 15b shows the estimated $C_a$ after 30 d of assimilation. All experiments approach the true value of $C_a$ (although starting from an overestimated value). When SIC and SIV are assimilated together, the IEnKF gives the best estimate of $C_a$, while it underestimates $C_a$ when only SIT and SIV are assimilated. The IEnKF also gives an increase in the estimate of $\alpha$ as shown in Fig. 15a. Although this demonstrates the difficulty in estimating both parameters, it is remarkable that it still leads to an improved forecast of model fields as shown in Fig. 16. This is known as a compensating effect (Bocquet, 2012). This apparent contradiction is related to the non-identifiability of the problem whereby the parameter estimation cannot fully reconstruct the true parameters. In practice, this implies that more than one set of parameters can produce the same observed fields. The identifiability issue can also be viewed from a physical standpoint: $C_a$ determines the motion of the sea ice via the wind forcing. With overestimated a priori $C_a$, the wind can cause fast sea-ice motions. $\alpha$ influences the rate at which sea ice transitions

from an elasto-brittle solid behaviour to a viscous fluid behaviour with an increasing level of damage. As discussed in Sect. 5.3, the increased $\alpha$ leads to more viscous sea ice that is more easily subjected to permanent deformations instead of fast transient elastic deformations. In this situation, the IEnKF controls the fast sea-ice motions by removing the elastic sea-ice deformations.

The identifiability issue can also be associated with incorrect cross-correlations between model fields and model parameters from different sources of errors as described in Eq. (7). We observe that assimilation correctly decreases $C_a$, showing a relatively accurate cross-correlation between the observations and $C_a$. This controls the impact of the external wind forcing, as a positively biased $C_a$ can amplify the uncertainty in the wind field, which increases the term $M^{u_a}\delta u_a$ in Eq. (7). As discussed in Sect. 5.3, the outstanding external wind forcing could lead to erroneous estimation of $\alpha$. Hence, after the external uncertainties are controlled, the ensemble should be able to develop reasonable cross-correlations between $\alpha$ and the observations. To test this, we adopt the following strategy in the 30 d assimilation experiments: (1) constraining the external uncertainty by estimating $C_a$ only for 10 d, (2) developing uncertainty from $\alpha$ by model forecasting without assimilation for 10 d and (3) estimating $\alpha$ for the last 10 d. This is labelled "seq" in Fig. 15. Although $C_a$ is still underestimated, we observe a reduction of error in $\alpha$. Also, as other experiments do not show a decreased $\alpha$ after an underestimated $C_a$, this may suggest the importance of step (2), where the uncertainty arising from $\alpha$ developed. We stress here that the 10 d is simply chosen for convenience and is not tuned.

To further consolidate our findings on the performance of the IEnKF in the simultaneous parameter estimations of $C_a$ and $\alpha$, extensive, and computationally demanding, investigations with different model states, forcings and parameters would be ideal. In fact, even though the model is a simplification of a full pan-Arctic sea-ice model, the computational power quickly scales up to beyond our resources.

With the aforementioned computational constraints in mind, we performed three additional experiments using different truths. In addition to the existing experiment (EXP$\alpha^L C_a^L$) where we assimilate all observations to estimate $\alpha$ and $C_a$, we implemented three other experiments where the truth uses (1) $\alpha = 4$ and $C_a = 3.5 \times 10^{-3}$ (EXP$\alpha^L C_a^H$), (2) $\alpha = 7$ and $C_a = 1.5 \times 10^{-3}$ (EXP$\alpha^H C_a^L$), and (3) $\alpha = 7$ and $C_a = 3.5 \times 10^{-3}$ (EXP$\alpha^H C_a^H$). Here, the superscript L and H denote that the truth is lower and higher respectively than the initial guess. We investigate relevant scenarios (e.g. under- or overestimation of the real values) in which the model's qualitative behaviour is of the same sort. The latter specifically concerns the model structural stability, that is to say the fact that the model is not subject to bifurcation of its general behaviour.

When the truth of model parameters is higher than the initial guess, the air drag coefficient is $2\sigma$ above the initial guess
<figure><img src="image_15.png" /><figcaption>**Figure 15.** The analysis of the $\alpha$ (a) and $C_a$ (b) parameter estimate based on a variety of combinations of observations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 16.** RMSE in the state estimate for scenario 4. (a) RMSE$^a$/RMSE$^{free}$. (b) RMSE$^f$/RMSE$^{free}$. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

while $\alpha$ is $1\sigma$ above the initial guess. Keeping the value up to $1\sigma$ higher than the initial guess of $\alpha$ is due to the changes in the dynamical regime of the model that occur beyond $1\sigma$. As discussed in Dansereau (2016), when $\alpha$ increases, the sea ice loses memory of the previous damage events, leading to increasing elasto-plastic behaviour. This implies that a different set of initial guesses might be needed for these dynamical regimes as we do not expect that the IEnKF to be able to provide reliable estimation when a change in dynamical regime occurs. Meanwhile, based on physical arguments, such a high $\alpha$ parameter is less likely to occur in reality where the true value is likely to be between 4 and 7. As shown in Fig. 17, $\text{EXP}\alpha^H C_a^L$ gives improvements in the air drag coefficient, and, though overestimated, the $\alpha$ parameter is increased as prescribed by the truth. In $\text{EXP}\alpha^L C_a^H$, improved $C_a$ estimation is obtained at the start of the estimation but deteriorates after 20 d. In $\text{EXP}\alpha^H C_a^H$, though the $C_a$ parameter is improved, the estimation of $\alpha$ is approximately 17, which is 1 order of magnitude larger than the truth. The deteriorated results occur when the true $C_a$ is higher than the initial guess, which corresponds to a strong wind forcing in the truth. One

possible explanation for this is that the correlations between $\alpha$ and the observed fields are not truthfully reflected due to the strong wind forcing. Nevertheless, if we first estimate the air drag coefficient, followed by a free forecast phase, and estimate $\alpha$ afterwards, the parameter estimation is improved. This shows that even if the same prior distribution is used, different dynamical regimes of the modelled truth can lead to different results in the experiments.

Our results suggest that the augmented state vector can have identifiability issues when both $C_a$ and $\alpha$ are biased. As shown in Table 5, the perturbation of $\alpha$ cannot provide greater forecast uncertainty compared to the positively biased $C_a$. We have proposed a strategy to overcome the issue with some discussions on alternative approaches. For example, a larger ensemble spread of $\alpha$ may increase the uncertainty in the forecast ensemble, which avoids a single dominant source of forecast uncertainty. We note that, based on our reasoning, the identifiability issue may not exist when $C_a$ is negatively biased. Nevertheless, this demonstrates a potential issue with the parameter estimation in MEB-type sea-ice models.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The analysis of the $\alpha$ (a) and $C_a$ (b) parameter estimate based on a variety of combinations of observations.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** RMSE in the state estimate for scenario 4. (a) RMSE$^a$/RMSE$^{free}$. (b) RMSE$^f$/RMSE$^{free}$. Columns display the individual physical model fields; rows refer to the (combination of) observations that are assimilated.</figcaption></figure>

while $\alpha$ is $1\sigma$ above the initial guess. Keeping the value up to $1\sigma$ higher than the initial guess of $\alpha$ is due to the changes in the dynamical regime of the model that occur beyond $1\sigma$. As discussed in Dansereau (2016), when $\alpha$ increases, the sea ice loses memory of the previous damage events, leading to increasing elasto-plastic behaviour. This implies that a different set of initial guesses might be needed for these dynamical regimes as we do not expect that the IEnKF to be able to provide reliable estimation when a change in dynamical regime occurs. Meanwhile, based on physical arguments, such a high $\alpha$ parameter is less likely to occur in reality where the true value is likely to be between 4 and 7. As shown in Fig. 17, $\text{EXP}\alpha^H C_a^L$ gives improvements in the air drag coefficient, and, though overestimated, the $\alpha$ parameter is increased as prescribed by the truth. In $\text{EXP}\alpha^L C_a^H$, improved $C_a$ estimation is obtained at the start of the estimation but deteriorates after 20 d. In $\text{EXP}\alpha^H C_a^H$, though the $C_a$ parameter is improved, the estimation of $\alpha$ is approximately 17, which is 1 order of magnitude larger than the truth. The deteriorated results occur when the true $C_a$ is higher than the initial guess, which corresponds to a strong wind forcing in the truth. One

possible explanation for this is that the correlations between $\alpha$ and the observed fields are not truthfully reflected due to the strong wind forcing. Nevertheless, if we first estimate the air drag coefficient, followed by a free forecast phase, and estimate $\alpha$ afterwards, the parameter estimation is improved. This shows that even if the same prior distribution is used, different dynamical regimes of the modelled truth can lead to different results in the experiments.

Our results suggest that the augmented state vector can have identifiability issues when both $C_a$ and $\alpha$ are biased. As shown in Table 5, the perturbation of $\alpha$ cannot provide greater forecast uncertainty compared to the positively biased $C_a$. We have proposed a strategy to overcome the issue with some discussions on alternative approaches. For example, a larger ensemble spread of $\alpha$ may increase the uncertainty in the forecast ensemble, which avoids a single dominant source of forecast uncertainty. We note that, based on our reasoning, the identifiability issue may not exist when $C_a$ is negatively biased. Nevertheless, this demonstrates a potential issue with the parameter estimation in MEB-type sea-ice models.
<figure><img src="image_17.png" /><figcaption>**Figure 17.** Simultaneous parameter estimation of $C_a$ and $\alpha$ with different sets of truths of these parameters. The horizontal red lines represent the truth values when they are either lower or higher than the initial guess. The inset in (a) is the estimated $\alpha$ in the experiment $\text{EXP}\alpha^H C_a^H$.</figcaption></figure>

## 5.5 Impact of parameter estimation on long-term prediction

It is worth studying the impact of the parameter estimation on the long-term performance of the model. To this end, we performed a 90 d long free run whereby the model parameters are “frozen” to the values obtained at the end of the DA period. The experiment is carried out for the case when both $C_a$ and $\alpha$ are estimated during the DA period using the configuration “seq” of Sect. 5.4. Results are displayed in Fig. 18.

Figure 18 shows that the RMSE of all observed model variables is smaller using the corrected model parameters compared to the free run that uses the initial guess. The improvement is persistent (in time) in SIV, while it increases for SIC and SIT due to the long timescale of these model variables. These results demonstrate the importance of correctly specified model parameters in long-term sea-ice forecasting.

## 6 Discussion

In this study, we investigate joint state and parameter estimation of an MEB sea-ice model using the IEnKF. This study focuses on the cross-covariance between model fields, which is crucial for correct parameter estimation. Given that the IEnKF is also a variational method, it minimises the cost

function in Eq. (1), and thus it adjusts the model parameter at each Gauss–Newton iteration. The adjusted, data-informed model parameters allow for a natural treatment of the nonlinearity thanks to the execution of the model within one analysis step or, had we used the IEnKF, within the entire assimilation window. The variational formulation also has other additional advantages over the traditional EnKF. For example, one can impose constraint optimisation and regularisation on the cost function to avoid numerical problems or to append physical constraints.

One important ingredient of the IEnKF is the use of an ensemble. The forecast ensemble of the IEnKF can suffer from ensemble collapse as discussed in Sect. 4.3, where inflation is used as an effective modification to the forecast ensemble. We adopt an adaptive inflation method, the IEnKF-N (Bocquet and Sakov, 2012). The adaptive inflation approach is shown to be applicable to complex models for joint state-parameter estimations with reasonable stopping criteria.

In addition to the constant ensemble spread, we also explored the use of an autoregressive model (NEA, 2007; Xie et al., 2017), instead of persistence, for the model parameters. However, the results (not shown) are elusive. We speculate that, as the autoregressive model resamples the ensemble uncertainty in the parameter at each time step, this introduces additional sampling error and the results may be subject to
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Simultaneous parameter estimation of $C_a$ and $\alpha$ with different sets of truths of these parameters. The horizontal red lines represent the truth values when they are either lower or higher than the initial guess. The inset in (a) is the estimated $\alpha$ in the experiment $\text{EXP}\alpha^H C_a^H$.</figcaption></figure>

## 5.5 Impact of parameter estimation on long-term prediction

It is worth studying the impact of the parameter estimation on the long-term performance of the model. To this end, we performed a 90 d long free run whereby the model parameters are “frozen” to the values obtained at the end of the DA period. The experiment is carried out for the case when both $C_a$ and $\alpha$ are estimated during the DA period using the configuration “seq” of Sect. 5.4. Results are displayed in Fig. 18.

Figure 18 shows that the RMSE of all observed model variables is smaller using the corrected model parameters compared to the free run that uses the initial guess. The improvement is persistent (in time) in SIV, while it increases for SIC and SIT due to the long timescale of these model variables. These results demonstrate the importance of correctly specified model parameters in long-term sea-ice forecasting.

## 6 Discussion

In this study, we investigate joint state and parameter estimation of an MEB sea-ice model using the IEnKF. This study focuses on the cross-covariance between model fields, which is crucial for correct parameter estimation. Given that the IEnKF is also a variational method, it minimises the cost

function in Eq. (1), and thus it adjusts the model parameter at each Gauss–Newton iteration. The adjusted, data-informed model parameters allow for a natural treatment of the nonlinearity thanks to the execution of the model within one analysis step or, had we used the IEnKF, within the entire assimilation window. The variational formulation also has other additional advantages over the traditional EnKF. For example, one can impose constraint optimisation and regularisation on the cost function to avoid numerical problems or to append physical constraints.

One important ingredient of the IEnKF is the use of an ensemble. The forecast ensemble of the IEnKF can suffer from ensemble collapse as discussed in Sect. 4.3, where inflation is used as an effective modification to the forecast ensemble. We adopt an adaptive inflation method, the IEnKF-N (Bocquet and Sakov, 2012). The adaptive inflation approach is shown to be applicable to complex models for joint state-parameter estimations with reasonable stopping criteria.

In addition to the constant ensemble spread, we also explored the use of an autoregressive model (NEA, 2007; Xie et al., 2017), instead of persistence, for the model parameters. However, the results (not shown) are elusive. We speculate that, as the autoregressive model resamples the ensemble uncertainty in the parameter at each time step, this introduces additional sampling error and the results may be subject to
<figure><img src="image_18.png" /><figcaption>**Figure 18.** RMSE and ensemble spread of the free run using prior parameters and parameters corrected by the DA when both $C_a$ and $\alpha$ are erroneous.</figcaption></figure>

the signal-to-noise ratio from the choice of the standard deviation of the autoregressive model. This may suggest the need for careful tuning for the autoregressive processes.

We also shed light on potential issues when the forecast uncertainty is driven mainly by the external wind field. In this case, the cross-covariance matrix reflects the error from the wind fields and the model fields. This is undesirable for the parameter estimations where we expect the error in the observed fields to be related to the parameter perturbations. This is particularly relevant for sea-ice DA where the wind field contributes to most of the uncertainties in the forecast ensemble. Meanwhile, we have also shown that the effect of external uncertainty depends on the model fields. The incorrect cross-correlations is more detrimental when the model fields are directly linked to the source of the uncertainty. For example, when only $\alpha$ is biased, assimilating only SIV shows severely problematic analysis. This may also suggest that coupled DA controlling the uncertainty in external forcing could improve the sea-ice parameter estimation.

## 7 Conclusions

We investigated the state and parameter estimation in a dynamics-only MEB sea-ice model under an idealised setup using the IEnKF (and the IEnKF-N). We mimicked the observation error and its spatial distribution with the forecast uncertainty driven primarily by the uncertainties in the wind field.

We adopted a fully multivariate approach whereby all model fields are estimated by DA utilising the cross-correlations between observations and model fields. Our results show that different combinations of the sea-ice observation fields can lead to different effects on the model state and parameter estimates. In general, it is useful to assimilate as many observations as possible. Potential issues with multivariate state and parameter estimation due to both the limitations of the DA method and the features of the sea-ice model dynamics are highlighted. We also demonstrate that these issues are surmountable and that successful multivariate state and parameter estimation using state-of-the-art ensemble DA approaches is possible.

Experiments in the perfect-model scenario show that, even if the sea-ice model is perfect, the limited ensemble size and uncertainties in the external forcing, e.g. the atmospheric wind, still limit the capability of DA to improve the sea-ice forecast especially for the unobserved model fields. We show that the forecast of SIV cannot be improved because it is strictly constrained by the wind field while other model fields with longer timescales show improved forecasts. This suggests that coupled DA that estimates the external forcing could improve the sea-ice forecast of model fields like SIV. In addition, the linear Gaussian assumption of DA methods can violate the bounds of various model fields in the sea-ice models, requiring post-processing of these fields. One potential treatment of the bound problem is anamorphosis, which can be applied analytically to the IEnKF as well as to other EnKF variants (e.g. Bocquet and Sakov, 2013; Simon and Bertino,
<figure><img src="image_21.png" /><figcaption>**Figure 18.** RMSE and ensemble spread of the free run using prior parameters and parameters corrected by the DA when both $C_a$ and $\alpha$ are erroneous.</figcaption></figure>

the signal-to-noise ratio from the choice of the standard deviation of the autoregressive model. This may suggest the need for careful tuning for the autoregressive processes.

We also shed light on potential issues when the forecast uncertainty is driven mainly by the external wind field. In this case, the cross-covariance matrix reflects the error from the wind fields and the model fields. This is undesirable for the parameter estimations where we expect the error in the observed fields to be related to the parameter perturbations. This is particularly relevant for sea-ice DA where the wind field contributes to most of the uncertainties in the forecast ensemble. Meanwhile, we have also shown that the effect of external uncertainty depends on the model fields. The incorrect cross-correlations is more detrimental when the model fields are directly linked to the source of the uncertainty. For example, when only $\alpha$ is biased, assimilating only SIV shows severely problematic analysis. This may also suggest that coupled DA controlling the uncertainty in external forcing could improve the sea-ice parameter estimation.

# 7 Conclusions

We investigated the state and parameter estimation in a dynamics-only MEB sea-ice model under an idealised setup using the IEnKF (and the IEnKF-N). We mimicked the observation error and its spatial distribution with the forecast uncertainty driven primarily by the uncertainties in the wind field.

We adopted a fully multivariate approach whereby all model fields are estimated by DA utilising the cross-correlations between observations and model fields. Our results show that different combinations of the sea-ice observation fields can lead to different effects on the model state and parameter estimates. In general, it is useful to assimilate as many observations as possible. Potential issues with multivariate state and parameter estimation due to both the limitations of the DA method and the features of the sea-ice model dynamics are highlighted. We also demonstrate that these issues are surmountable and that successful multivariate state and parameter estimation using state-of-the-art ensemble DA approaches is possible.

Experiments in the perfect-model scenario show that, even if the sea-ice model is perfect, the limited ensemble size and uncertainties in the external forcing, e.g. the atmospheric wind, still limit the capability of DA to improve the sea-ice forecast especially for the unobserved model fields. We show that the forecast of SIV cannot be improved because it is strictly constrained by the wind field while other model fields with longer timescales show improved forecasts. This suggests that coupled DA that estimates the external forcing could improve the sea-ice forecast of model fields like SIV. In addition, the linear Gaussian assumption of DA methods can violate the bounds of various model fields in the sea-ice models, requiring post-processing of these fields. One potential treatment of the bound problem is anamorphosis, which can be applied analytically to the IEnKF as well as to other EnKF variants (e.g. Bocquet and Sakov, 2013; Simon and Bertino,
2009). We did not explore this avenue here, but it may constitute an interesting follow-up work; note however that the majority of sea-ice DA does not use anamorphosis, since it can cause numerical imbalances (Bocquet and Sakov, 2013). Nevertheless, with suitable localisation and sufficient observations, we show improvements for all model fields, both observed and unobserved fields (e.g. stress, cohesion and damage).

We choose two global parameters to be estimated, the air drag coefficient, $C_a$, and the damage parameter, $\alpha$. The air drag coefficient is closely related to the external uncertainties from the wind field, while $\alpha$ affects the mechanical and dynamical regime of the sea ice. In this study, we choose one specific setup of the model parameters as the truth based on experiments in Dansereau (2016). With our chosen model setup, in a model with only one biased parameter, the DA can reduce the parameter bias and improve the model forecast. However, in the case where only SIV is assimilated, our results show that it can lead to model imbalance. Also, as it is closely related to the external uncertainties, the cross-correlation between SIV and $\alpha$ can be incorrectly specified. This shows the importance of assimilating multiple sea-ice observations and the potential difficulty in sea-ice DA when the ensemble uncertainties are primarily driven by the external uncertainties (wind field).

When both model parameters are biased, an identifiability problem arises. This highlights the caveat that the ensemble spread can come from different sources of uncertainties rather than purely from uncertainties in the parameters. When one source of uncertainty dominates the ensemble uncertainty, the estimated model parameters can deviate from the truth even if a better forecast is achieved for those observed fields. We proposed a strategy that can mitigate such an issue in our specific test case by first controlling the dominating error from the external forcing and later estimating the sea-ice internal parameters.

A number of open questions are still at stake. For example, the MEB model contains multiple dynamical regimes as the sea ice can be viewed as different materials based on the state of sea ice. It is of theoretical interest to investigate the parameter estimations under different dynamical regimes of the MEB model. We observed improvements in the fully multivariate update, but it is still unclear whether these improvements can be observed in the full Arctic sea-ice predictions compared to a state-of-the-art operational setup. Another point that should be addressed in future work is the development of a rigorous approach to handle bounded variables like concentration and damage.

On the side of the ensemble generation, in this study, the wind forcing serves as the sole source of uncertainty. Nevertheless, in operational sea-ice models, multiple potential uncertain external forcing sources are present. A study of MEB-like models’ sensitivity to different external forcing appears to be another avenue worth pursuing, along the lines of the

analysis of sensitivity to external wind and cohesion parameters in Cheng et al. (2020).

We finally mention the current development of data assimilation approaches specifically adapted to models using a discontinuous Galerkin dynamical core, as in the sea-ice model used in this study (Pasmans et al., 2023). These novel methods may pave the way to solving some of the aforementioned issues, and the authors are currently investigating their use in the present context.

## Appendix A: Model parameters

Here we present the model parameters used in the modelled truth in Table A1.

<table><caption>Table A1. The default model parameters used in the modelled truth.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Value</th></tr></thead><tbody><tr><td>Poisson's ratio</td><td>ν</td><td>0.3</td></tr><tr><td>Internal friction coefficient</td><td>μ</td><td>0.7</td></tr><tr><td>Ice density</td><td>ρ</td><td>900 kg m<sup>-3</sup></td></tr><tr><td>Elastic (shear) wave propagation speed</td><td>c</td><td>500 m s<sup>-1</sup></td></tr><tr><td>Damage parameter</td><td>α</td><td>4.0</td></tr><tr><td>Undamaged elastic modulus</td><td>E<sup>0</sup></td><td>2c<sup>2</sup>(1 + ν)ρ</td></tr><tr><td>Undamaged relaxation time (1/τ<sup>0</sup>/E<sup>0</sup>)</td><td>λ<sup>0</sup></td><td>10<sup>7</sup> s</td></tr><tr><td>Undamaged apparent viscosity</td><td>η<sup>0</sup></td><td>λ<sup>0</sup>E<sup>0</sup></td></tr><tr><td>Minimum apparent viscosity</td><td>η<sub>min</sub></td><td>10<sup>4</sup> Pa s</td></tr><tr><td>Minimum cohesion</td><td>C<sub>min</sub></td><td>5000 Pa</td></tr><tr><td>Model time step</td><td>Δt</td><td>30 s</td></tr><tr><td>Mean model resolution</td><td>Δx</td><td>15 km</td></tr><tr><td>Characteristic time for damage</td><td>t<sub>d</sub></td><td>Δt</td></tr><tr><td>Characteristic time for healing</td><td>t<sub>h</sub></td><td>5 × 10<sup>5</sup> s</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub></td><td>1.3 kg m<sup>-3</sup></td></tr><tr><td>Air drag coefficient</td><td>C<sub>d<sub>a</sub></sub></td><td>1.5 × 10<sup>-3</sup></td></tr><tr><td>Water density</td><td>ρ<sub>w</sub></td><td>1027 kg m<sup>-3</sup></td></tr><tr><td>Water drag coefficient</td><td>C<sub>d<sub>w</sub></sub></td><td>5.5 × 10<sup>-3</sup></td></tr><tr><td>Parameter used in coupling E and η to A</td><td>c*</td><td>20</td></tr></tbody></table>

## Appendix B: Wind field

The wind field is prescribed as a series of passing cyclonic storms combined with a constant background wind. To have control over the vorticity and divergence of the wind field, we use the Helmholtz decomposition:

$$ \mathbf{u}_a = \nabla \Phi + \mathbf{k} \times \nabla \Psi, \qquad (B1) $$

where $\mathbf{u}_a = (u_a, v_a)$, $\Phi$ is the velocity potential, and $\Psi$ is the streamfunction with the solenoidal wind component $\mathbf{u}_a^s = \nabla \Phi$ and the divergent wind component $\mathbf{u}_a^d = \mathbf{k} \times \nabla \Psi$. The
2009). We did not explore this avenue here, but it may constitute an interesting follow-up work; note however that the majority of sea-ice DA does not use anamorphosis, since it can cause numerical imbalances (Bocquet and Sakov, 2013). Nevertheless, with suitable localisation and sufficient observations, we show improvements for all model fields, both observed and unobserved fields (e.g. stress, cohesion and damage).

We choose two global parameters to be estimated, the air drag coefficient, $C_a$, and the damage parameter, $\alpha$. The air drag coefficient is closely related to the external uncertainties from the wind field, while $\alpha$ affects the mechanical and dynamical regime of the sea ice. In this study, we choose one specific setup of the model parameters as the truth based on experiments in Dansereau (2016). With our chosen model setup, in a model with only one biased parameter, the DA can reduce the parameter bias and improve the model forecast. However, in the case where only SIV is assimilated, our results show that it can lead to model imbalance. Also, as it is closely related to the external uncertainties, the cross-correlation between SIV and $\alpha$ can be incorrectly specified. This shows the importance of assimilating multiple sea-ice observations and the potential difficulty in sea-ice DA when the ensemble uncertainties are primarily driven by the external uncertainties (wind field).

When both model parameters are biased, an identifiability problem arises. This highlights the caveat that the ensemble spread can come from different sources of uncertainties rather than purely from uncertainties in the parameters. When one source of uncertainty dominates the ensemble uncertainty, the estimated model parameters can deviate from the truth even if a better forecast is achieved for those observed fields. We proposed a strategy that can mitigate such an issue in our specific test case by first controlling the dominating error from the external forcing and later estimating the sea-ice internal parameters.

A number of open questions are still at stake. For example, the MEB model contains multiple dynamical regimes as the sea ice can be viewed as different materials based on the state of sea ice. It is of theoretical interest to investigate the parameter estimations under different dynamical regimes of the MEB model. We observed improvements in the fully multivariate update, but it is still unclear whether these improvements can be observed in the full Arctic sea-ice predictions compared to a state-of-the-art operational setup. Another point that should be addressed in future work is the development of a rigorous approach to handle bounded variables like concentration and damage.

On the side of the ensemble generation, in this study, the wind forcing serves as the sole source of uncertainty. Nevertheless, in operational sea-ice models, multiple potential uncertain external forcing sources are present. A study of MEB-like models' sensitivity to different external forcing appears to be another avenue worth pursuing, along the lines of the

analysis of sensitivity to external wind and cohesion parameters in Cheng et al. (2020).

We finally mention the current development of data assimilation approaches specifically adapted to models using a discontinuous Galerkin dynamical core, as in the sea-ice model used in this study (Pasmans et al., 2023). These novel methods may pave the way to solving some of the aforementioned issues, and the authors are currently investigating their use in the present context.

## Appendix A: Model parameters

Here we present the model parameters used in the modelled truth in Table A1.

<table><caption>Table A1. The default model parameters used in the modelled truth.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Value</th></tr></thead><tbody><tr><td>Poisson's ratio</td><td>ν</td><td>0.3</td></tr><tr><td>Internal friction coefficient</td><td>μ</td><td>0.7</td></tr><tr><td>Ice density</td><td>ρ</td><td>900 kg m<sup>-3</sup></td></tr><tr><td>Elastic (shear) wave propagation speed</td><td>c</td><td>500 m s<sup>-1</sup></td></tr><tr><td>Damage parameter</td><td>α</td><td>4.0</td></tr><tr><td>Undamaged elastic modulus</td><td>E<sup>0</sup></td><td>2c<sup>2</sup>(1 + ν)ρ</td></tr><tr><td>Undamaged relaxation time (1/τ<sup>0</sup>/E<sup>0</sup>)</td><td>λ<sup>0</sup></td><td>10<sup>7</sup> s</td></tr><tr><td>Undamaged apparent viscosity</td><td>η<sup>0</sup></td><td>λ<sup>0</sup>E<sup>0</sup></td></tr><tr><td>Minimum apparent viscosity</td><td>η<sub>min</sub></td><td>10<sup>4</sup> Pa s</td></tr><tr><td>Minimum cohesion</td><td>C<sub>min</sub></td><td>5000 Pa</td></tr><tr><td>Model time step</td><td>Δt</td><td>30 s</td></tr><tr><td>Mean model resolution</td><td>Δx</td><td>15 km</td></tr><tr><td>Characteristic time for damage</td><td>t<sub>d</sub></td><td>Δt</td></tr><tr><td>Characteristic time for healing</td><td>t<sub>h</sub></td><td>5 × 10<sup>5</sup> s</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub></td><td>1.3 kg m<sup>-3</sup></td></tr><tr><td>Air drag coefficient</td><td>C<sub>d<sub>a</sub></sub></td><td>1.5 × 10<sup>-3</sup></td></tr><tr><td>Water density</td><td>ρ<sub>w</sub></td><td>1027 kg m<sup>-3</sup></td></tr><tr><td>Water drag coefficient</td><td>C<sub>d<sub>w</sub></sub></td><td>5.5 × 10<sup>-3</sup></td></tr><tr><td>Parameter used in coupling E and η to A</td><td>c*</td><td>20</td></tr></tbody></table>

## Appendix B: Wind field

The wind field is prescribed as a series of passing cyclonic storms combined with a constant background wind. To have control over the vorticity and divergence of the wind field, we use the Helmholtz decomposition:

$$ \mathbf{u}_a = \nabla \Phi + \mathbf{k} \times \nabla \Psi, \qquad (B1) $$

where $\mathbf{u}_a = (u_a, v_a)$, $\Phi$ is the velocity potential, and $\Psi$ is the streamfunction with the solenoidal wind component $\mathbf{u}_a^s = \nabla \Phi$ and the divergent wind component $\mathbf{u}_a^d = \mathbf{k} \times \nabla \Psi$. The