Articles | Volume 14, issue 12
Research article
04 Dec 2020
Research article |  | 04 Dec 2020

Parameter optimization in sea ice models with elastic–viscoplastic rheology

Gleb Panteleev, Max Yaremchuk, Jacob N. Stroh, Oceana P. Francis, and Richard Allard

The modern sea ice models include multiple parameters which strongly affect model solution. As an example, in the CICE6 community model, rheology and landfast grounding/arching effects are simulated by functions of the sea ice thickness and concentration with a set of fixed parameters empirically adjusted to optimize the model performance. In this study, we consider the extension of a two-dimensional elastic–viscoplastic (EVP) sea ice model using a spatially variable representation of these parameters. The feasibility of optimization of the landfast sea ice parameters and rheological parameters is assessed via idealized variational data assimilation experiments with synthetic observations of ice concentration, thickness and velocity. The experiments are configured for a 3 d data assimilation window in a rectangular basin with variable wind forcing.

The tangent linear and adjoint models featuring EVP rheology are found to be unstable but can be stabilized by adding a Newtonian damping term into the adjoint equations. A set of observation system simulation experiments shows that landfast parameter distributions can be reconstructed after 5–10 iterations of the minimization procedure. Optimization of sea ice initial conditions and spatially varying parameters in the stress tensor equation requires more computation but provides a better hindcast of the sea ice state and the internal stress tensor. Analysis of inaccuracy in the wind forcing and errors in sea ice thickness observations show reasonable robustness of the variational DA approach and the feasibility of its application to available and incoming observations.

1 Introduction

Due to the significant decline of sea ice volume and the increase in maritime activity in the Arctic Ocean over the past decades, an accurate sea ice hindcast/forecast has become an important component of the data assimilation systems for the region. Currently, there are several community sea ice models broadly used for modeling and/or reconstruction of the Arctic Ocean state through various data assimilation (DA) algorithms. Many of these models (e.g., Heimbach et al., 2010; Zhang and Rothrock, 2003; Vancoppenolle et al., 2009; Massonnet et al., 2015) are based on the viscoplastic (VP) rheology proposed by Hibler (1979). In the last decades, several numerical approaches have been proposed for solving the VP problem (Hibler, 1979; Zhang and Hibler, 1997; Lemieux et al., 2008). These approaches are based on implicit solvers and require a significant number of iterations to achieve full convergence. Application of more efficient quasi-Newtonian solvers suffers from the lack of robustness, and these are usually applied in the sea ice models of intermediate resolution (Lemieux et al., 2012; Losch et al., 2014). In addition, implicit VP solvers are less convenient for implementing on massively parallel supercomputer architectures. Despite these inconveniences, the currently known four-dimensional variational (4D-Var) sea ice models employ an implicit VP solver for time integration of the tangent linear and adjoint (TLA) models (Kauker et al., 2009; Heimbach et al., 2010).

The elastic–VP (EVP) rheology was proposed as an alternative explicit method which can be easily adopted for supercomputer architectures (Hunke and Duckowicz, 1997). In addition to the VP parameterization of the internal stress tensor, EVP includes an additional elastic term which requires internal subcycling to damp elastic oscillations in order to achieve the VP solution. The most popular sea ice model with EVP rheology is the Community sea ICE model (CICE, Hunke et al., 2010), which is currently maintained and developed by a group of institutions in North America and Europe known as the CICE Consortium. This model is widely used in sea ice and coupled sea–ice modeling (e.g., Posey et al., 2010, Metzger et al., 2014, Yaremchuk et al., 2019), and there are multiple examples of 3D-Var, ensemble-based DA systems utilizing the CICE model (e.g., Zhang and Bitz, 2018). However, to our knowledge, there is no 4D-Var DA system based on the CICE model yet.

CICE6 includes several parameterizations of the sea ice rheology including the formulation of Hibler (1979). This parameterization includes three major parameters (P*, e and α), describing, respectively, the dimensional maximum ice strength per unit thickness, the ratio of yield ellipse major axes and the nondimensional scaling of ice strength with its compactness. While this parameterization is not the default option in CICE5/6, it is still widely used in sea ice modeling and DA applications at the Naval Research Laboratory.

For modeling landfast ice near the coast and in straits (the location of a so-called arching phenomena), an additional parameter kT has been introduced to model the tensile strength (König Beatty and Holland, 2010). This parameter is absent in the traditional (i.e., Hibler, 1979) elliptical yield curve formulation. Lemieux et al. (2015) proposed a number of additional parameters (k1,k2 and αb) for better parameterization of the landfast ice grounded in the shallow regions. Formally, these landfast ice parameters are not related to the sea ice rheology. To simplify the presentation, we shall refer the entire set {P*,e,kt,k1,k2} as rheological parameters (RPs), suggesting that all of them influence sea ice dynamics.

Typically, the abovementioned rheological parameters are constants, and their values are defined empirically from multiple numerical experiments. RPs such as P* and e reflect the model parameterization rather than physics and are not directly observable (Kreyscher et al., 2000) but are nevertheless known to range within certain limits (Harder and Fischer 1999). As a few examples, the typical values of P* determined from sea ice drift were diagnosed to vary within 27.5 kN/m2 (Hibler and Walsh, 1982), 15–20 kN/m2 (Kreyscher et al., 1997, 2000) and 30–45 kN/m2 (Tremblay and Hakakian, 2006). These studies indicate the existence of significant variations of P* estimates, which may be attributed to both nonphysical considerations (such as spatially variable model resolution) and spatiotemporal variations of Arctic sea ice.

The numerical experiments of Lemieux et al. (2016) using a coarse-resolution pan-Arctic CICE-NEMO model have shown that kT=0.2 provides the best agreement with landfast strength observations in the Kara Sea, when the ellipse axes ratio ranges within 1.2–1.4. The most sensitive parameters for sea ice grounding are the critical thickness parameter k1 and maximum basal stress parameter k2. The optimal values were found to be 8 and 15 N/m3, respectively, with higher sensitivity with respect to k1 (Lemieux et al., 2015). However, fixed values of kT,k1 and k2 cannot provide a universally good performance of landfast modeling for different parts of the Arctic Ocean, suggesting that these parameters are a function of local environmental conditions. The physical properties of sea ice also strongly depend on sea ice salinity, temperature and ice age (e.g., Anderson and Weeks, 1958), indicating that rheological properties may vary between different sea ice categories.

Thus, numerous modeling experiments and sea ice observations (e.g., Juricke et al., 2013; Toyota and Kimura, 2018) indicate that spatially varying and properly optimized RPs should significantly improve the sea ice model performance. There were multiple attempts to define sea ice model parameters in an optimal way. The early attempts followed a traditional trial-and-error approach, involving multiple runs of a sea ice model with different sets of RPs (e.g., Miller et al., 2006; Uotila et al., 2012), while others utilized more advanced methods based on the Green's function approach (Nguyen et al., 2011), ensemble Kalman filtering (Massonnet et al., 2014) and genetic algorithms (Sumata et al., 2019). However, all of these attempts sought a relatively small set of spatially invariant sea ice model parameters in order to provide an optimal sea ice model solution for a period of several years or decades. The application of these algorithms for optimizing spatially varying RPs was not considered and, from our point of view, is not straightforward due to high computational overhead. Also note that the above-listed algorithms are not suitable for simultaneous optimization of other model parameters such as initial and open-boundary conditions, or external forcing fields.

The major objective of our study is to find a numerically feasible method for simultaneous optimization of multiple model parameters including the spatially varying RPs, initial conditions and forcing fields in the sea ice models based on EVP solvers (e.g., CICE model). The framework of a Naval Research Laboratory research project to identify spatially varying landfast ice parameters in the CICE6 model guided the priorities and objectives of this study. We suggest that, if successful, this approach can be adopted to optimize RPs in operational sea ice models (e.g., Posey et al., 2010; Metzger et al., 2014) and provide a more accurate short-term (3–7 d) sea ice forecast.

The feasibility of reconstructing spatially varying fields of P* and e (as well as other model parameters) through the variational assimilation of synthetic observations of sea ice velocity/concentration/thickness (SIV/SIC/SIT) was recently analyzed by Stroh et al. (2019) in the framework of a one-dimensional (1D) (zonal) VP sea ice model. It was found that variational DA allows for a reasonable reconstruction of spatially varying P* and e in regions with strong convergence and significantly improves short-range hindcast/forecast of the sea ice state. In particular, it was shown that optimization of spatially varying P* and e provides more accurate reconstruction of ridging areas, which cannot be achieved by optimizing the initial sea ice state only.

Note that optimization of RPs through the 4D-Var DA approach allows us to efficiently use all available sea ice observations, including observations of sea ice velocity, which are rarely used for assimilation in sea ice DA systems controlled by the initial conditions only. This is due to weak sensitivity of the sea ice state with respect to sea ice velocity initial conditions (e.g., Kauker, et al., 2009). Augmenting the 4D-Var control vector with RPs allows us to use sea ice velocity observations for consistent adjustment of the RPs and/or atmospheric forcing, providing a better sea ice forecast (Stroh et al., 2019).

In this study, we extend the investigation to analyze the feasibility of RP optimization within a more advanced two-dimensional (2D) sea ice model based on the EVP rheology formulation of Lemieux et al. (2016). Our analysis is based on application of the 4D-Var and Observing System Simulation Experiment (OSSE) (e.g., Nitta, 1969; Arnold and Dey, 1986; Nichols 2003, 2010) and follows the conventional twin-data experiment approach (e.g., Goldberg and Heimbach, 2013).

Similarly to Stroh et al. (2019), we developed the corresponding EVP TLA models and analyzed the feasibility of optimizing spatially varying ellipse ratio and sea ice compressive strength. In addition, we also analyzed the effects of optimizing two of the landfast sea ice parameters introduced by Lemieux et al. (2016). Through multiple OSSEs, we evaluate the quality of the RP reconstruction and analyze the impact of spatially varying RPs on the sea ice state. A similar approach was recently proposed for the optimization of the basal stress parameters in an ice sheet model (Goldberg and Heimbach, 2013).

Currently, satellite sea ice observations are typically available daily with a reasonably dense spatial resolution. Analysis of synthetic-aperture radar (SAR) images (e.g., Panteleev et al., 2019) indicates that in the marginal sea ice zone, pancake/cake ice with floe sizes of  1–20 m may be easily replaced by floes exceeding 1 km in size in 1 week. As a consequence, we configured the OSSEs with a 3 d DA window assuming that sea ice features do not move very far from their initial position during this period. This assumption suggests that such an approach should have more impact on the short-term sea ice forecast.

The paper is organized as follows: Sect. 2 describes the implemented sea ice model, the details of the TLA codes, and the generation of synthetic observations and the first-guess solution used in OSSEs. Results of these experiments are described in Sects. 3 (optimization of landfast sea ice parameters) and 4 (optimization of compressive strength and ellipse ratio), with a special focus on the feasibility of optimizing spatially varying RPs in the context of present and future observational coverage of sea ice at high latitudes. Section 5 summarizes the work and discusses directions of future research.

2 Sea ice model and its 4D-Var implementation

This section provides details of the sea ice model formulation and its associated linearizations, outlines the variational data assimilation system used for optimizing model parameters, and describes synthetic observations used to do so. To distinguish between the parameter fields spatially varying in 2D and the fixed parameter values that were not subject to optimization, the latter are marked by tildes. The major parameters of the model are listed in Table 1.

2.1 EVP sea ice model

2.1.1 Formulation

In the present study, we employed the sea ice model formulation of Lemieux et al. (2016) with the basal stress parameterization and generalized Hibler (1979) yield curve. Equations of the model describe EVP ice physics coupled with sea ice dynamics, which is forced by the stresses τ exerted on ice through its interaction with the bottom τb, atmosphere τa and the ocean τw:


Here div and tr are the divergence and trace operators; k is the vertical unit vector; I is the 2 × 2 identity matrix; h,A and u={u,v} are the 2D fields of ice effective thickness, concentration and velocity; and σ and ϵ˙ are the 2D fields of ice stress and the deformation rate tensors:

(5) σ = σ x x σ x y σ x y σ y y ; ϵ ˙ = 1 2 2 x u x v + y u x v + y u 2 y v .

The scalar field Δ used for normalizing the right-hand side (rhs) in (2) is computed using the following expression (e.g., Hunke, 2001):

(6) Δ ( ϵ ˙ ) = 1 e ( e 2 - 1 ) ( tr ϵ ˙ ) 2 + 2 tr ( ϵ ˙ 2 ) 1 / 2 .

To avoid numerical singularities at ϵ˙=0, the values of Δ are limited from below by the additional parameter Δ̃=10-10 s−1, so that Δ(ϵ˙)=max(Δ,Δ̃).

The empirical parameters T̃d,P*,kT and e in Eqs. (1)–(4) define the elastic damping scale of sea ice, its internal pressure, isotropic tensile strength and the yield curve axes ratio, respectively.

The internal ice pressure Pp is related to ice thickness and concentration in accordance with the rheology of Hibler (1979):

(7) P p = P * h A exp - α ̃ ( 1 - A ) .

The typical values of the ice strength parameter P* and α̃ are listed in Table 1.

The bottom and ocean stresses in Eq. (1) were parameterized in accordance with Lemieux et al. (2015, 2016) and Hunke and Lipscomb (2008):


where θ is the Heaviside step function, hb is the ocean depth, RΘ is the 2×2 matrix rotating the velocity vector by the turning angle Θ counterclockwise, ũ0=10-5 m/s and uw is the water velocity (set to zero in the present study). The values of other parameters (C̃w,ρw,α̃b,k̃1 and k2) are listed in Table 1.

Table 1Model and assimilation system configuration parameters.

Download Print Version | Download XLSX

In contrast to the previous studies, where the free empirical parameters P*,e,kT and k2 were assumed to be constant, the present study attempts to retrieve their spatial variability from synthetic (satellite) observations of the sea ice state vector C{h,A,u} using the variational DA technique.

2.1.2 Numerical scheme

Numerical formulation of the model closely follows the EVP numerics given in Lemieux et al. (2016).

Introducing notations σ1=σxx+σyy, σ2=σxx-σyy, σ3=σxy, ϵ˙1=xu+yv, ϵ˙2=xu-yv, ϵ˙3=yu+xv, bulk viscosity ζ=Pp(1+kT)/2Δ, and replacement pressure P=PpΔ/Δ, Eq. (2) can be split into three decoupled relationships:


The first equation is obtained by taking the trace of Eq. (2). Subtracting the equation for σyy from the one for σxx in the set of four relationships (2) yields Eq. (11), while Eq. (12) is just the equation for σxy extracted from the same set. Equations (10)–(12) are advanced in time s using the Euler scheme with the subcycling time step δts=εT̃d (ε=0.00347; see Table 1):


Hereinafter, all the fields without superscripts are taken from the previous time steps. After that, velocity components in Eq. (1) are advanced in time:


Here m=ρ̃hA and τa is the atmospheric forcing. After the Eqs. (13)–(17) are advanced in time for 400 time steps, the ice thickness and concentration fields are updated using the obtained velocity un and a simplified Lax–Wendroff scheme:


where D^ is the Laplacian operator and ν=δtun2/2 is the stabilizing diffusion coefficient. All spatial derivatives present in Eqs. (13)–(19) were discretized by finite differences on the Arakawa B-grid. At the rigid boundaries, we used the condition u=0. Initial conditions for u were set either to zero or defined through the model integration for 1 h. Initial conditions for A and h were specified as arbitrary functions.

The two-stage time stepping EVP procedure, described above, differs from the VP formulation in that iterations of the implicit solver of the VP formulation are replaced by the explicit adjustment of the stress tensor components (Eqs. 1315) in the internal time loop.

2.2 Variational DA with EVP sea ice model

2.2.1 Strong constraint formulation

In the variational DA experiments we used the so-called strong constraint state-space formulation of the problem, which minimizes a user-defined cost function on the manifold whose structure is specified by the equations of the model. The cost function 𝒥 was defined by

(20) J = 1 2 Ω W h ( h - h ) 2 + W A ( A - A ) 2 + W u ( u - u ) 2 + W ̃ h ( D ^ h ) 2 + W ̃ A ( D ^ A ) 2 + W ̃ u ( D ^ u ) 2 .

Here, W and W̃ denote the nonzero elements of the user-defined (diagonal) inverse error covariance matrices of the fields in the round brackets, simulated observations are denoted by primes, and summation is made over the entire space–time computational grid Ω. Note that the first three terms attract the optimized solution to the data, while the last three tend to penalize grid-scale components and enforce smoothness of the optimized fields.

In order to constrain the minimization process to the manifold M defined by the model Eqs. (13)–(19), we introduce notation X for a vector of state variables in all the grid points of Ω and define the vector of control variables C=[C0,Cp], which includes the initial state of the model C0X|t=0, and other control fields Cp, which contain rheological parameters, and atmospheric forcing errors. Note that the vector of model trajectory X is a nonlinear function of the control vector C, whose constituent Cp was defined on a sparser 2D grid (in every 15th node of the computational grid for most of the conducted experiments) and then spatially interpolated on the model grid using the bilinear interpolation operator .

To constrain the minimization to M(X,C), we introduce the vector of Lagrangian multipliers X^ (adjoints of the state variables X) in Ω (e.g., Le Dimet and Talagrand, 1986) and minimize the modified cost function

(21) J m = J + J ^ J ( X ) + M ( X , C ) T X ^ ,

with respect to X,X^ and C, where 𝖳 denotes the transposition. The minimum point is defined by setting the gradients of 𝒥m (i.e., variations of 𝒥m that are linear in δX,δX^ and δC) to zero:


As it is seen, Eq. (22) simply presents the nonlinear model trajectory specified by a given set of control variables. The second relationship contains the derivatives of 𝒥 and J^ with respect to X and involves the product of the model operator MX linearized in the vicinity of X (the tangent space to M at X, or the TL model) by the vector of respective perturbations δX. As soon as the current iteration of X is determined by running the nonlinear model, Eq. (23) can be solved by running the adjoint model (transpose of MX) forced by the derivatives of 𝒥 to obtain the values of the adjoint variables X^. Finally, the derivatives of 𝒥m are computed from Eq. (24) using the chain rule and the values of X and X^ derived from solving Eqs. (22)–(23).

The minimization of the cost function with respect to C was performed using the limited-memory quasi-Newtonian (M1QN3) method of Gilbert and Lemaréchal (1989) with the additional range constraints for the selected control fields (Sect. 2.3) performed after each iteration. The M1QN3 algorithm employs the above procedure for computation of the cost function gradient δ𝒥∕δC with respect to C for a given value of the control fields.

To sum up, the minimization procedure can be outlined as follows:

  1. specify a control vector C={u(0,x),A(0,x),h(0,x),Cp(t,x)} at the nth iteration;

  2. run the forward model (Eqs. 1319) to compute Xn, the derivatives δ𝒥∕δX and the value of 𝒥n (Eq. 20);

  3. run the adjoint model forced by -δJ/δX (cf. Eq. 23) to obtain X^n;

  4. compute 𝒥mδC (Eq. 24);

  5. update the control variables using the M1QN3 software;

  6. repeat steps (1)–(5) until convergence of the M1QN3 minimization procedure.

Technically, apart from developing the model code (Eqs. 1319), the outlined optimization requires development of the routines for computing the gradients as well as the tangent linear model MX and its adjoint MXT. The machinery of deriving these codes is based on the rules of differentiation and was realized in multiple software packages (e.g., Giering and Kaminski, 1998; AutoDiff:, last access: 1 December 2020; OpenAD:, last access: 1 December 2020; Goldberg and Heimbach, 2013).

More details on the variational techniques of data assimilation in different geophysical applications can be found in numerous publications (e.g., Penenko, 1981; Le Dimet, 1982; Lewis and Derber, 1985; Le Dimet and Talagrand, 1986; Wunsch, 1996; Errico, 1997). Note that both finite-difference TL and adjoint models are completely defined by the finite-difference scheme of the forward model, thus allowing application of the abovementioned (semi-)automatic TL/adjoint compilers. An alternative is used in our implementation: the code for MX is derived by an analytical linearization of the discretized forward model, while the adjoint code (Eq. 23) is obtained by analytical differentiation of 𝒥m with respect to the argument of the TL code. However, numerical stability of the nonlinear forward model does not guarantee stability of the respective TL and adjoint models, and it requires proper regularization (e.g., Hoteit et al., 2005) to move the eigenvalues of unstable eigenmodes of MX inside the unit circle. This numerical issue is addressed in the following section.

2.2.2 Adjoint and tangent linear models

The TL code was derived by analytic differentiation of the abovementioned numerical scheme in the vicinity of the finite-difference model trajectory. The adjoint code was obtained by implicit transposition of the sparse matrix in the code simulating the action of the TL operator MX on a perturbed state vector δX. Similar to the TL derivation, this procedure was performed by analytic differentiation with respect to δX of the code for computing X^TMXδX (cf. Eq. 23). More detailed description of the TLA codes and the gradients with respect to the control variables can be found in Appendix A.

The most laborious part of deriving the TLA codes was associated with linearizing the rhs of Eq. (2) with respect to ice velocities and RPs. Note that the first term in the rhs of the linearized Eq. (2) is proportional to the first derivatives of the velocity perturbations δu. As a consequence, the components of σ are linear in the first derivatives of δu after taking the explicit time step δts in the linearized Eq. (2). Moreover, the first-order derivatives in u keep their presence in the rhs of the linearized Eq. (1) due to spatial variability of the background fields in Eqs. (6) and (7).

This property of the TL equations of the subsystem (1)–(2) may require additional care when specifying the subcycling time step δts because the gradients of the background fields of h and A may invoke considerably larger propagation speeds of the effective elastic waves than those present in the original nonlinear model. Consequently, the TL code could be constrained by a more stringent stability criterion and require even smaller subcycling time steps than those used in the integration of the full nonlinear model. In particular, the nonlinear stability criterion could be violated, for example, in areas of strong ice convergence. In such regions, Eq. (7) implies that large SIT gradients may generate large coefficients of first-order derivatives of δu in the TL code for the second term in the rhs of Eq. (2).

Preliminary numerical experimentation with the TL code exposed a necessity to reduce δts as the TL solutions demonstrated uncontrollable amplification of velocity perturbations over the areas of strong sea ice convergence in the background fields. Our attempts to reduce δts by an order in magnitude reduced this type of amplification with a limited success. A similar instability of the TL EVP solver has been observed in the MITgcm sea ice model (Martin Losch, personal communication, 2019).

Instability of the linearized codes in the strongly nonlinear regimes of the parent model is a well-known phenomenon in the ocean general circulation models (OGCMs). A heuristic solution to this problem was proposed by Hoteit et al. (2005), who added an extra diffusion in TLA codes to suppress unstable small-scale harmonics. This kind of treatment is achieved, however, at the expense of reducing the TLA code accuracy (e.g., Yaremchuk et al., 2009). Later, Yaremchuk and Martin (2014) established a connection between the length of the DA window and the magnitude of the diffusion tensor in the TLA regularization terms.

However, in the sea ice model considered, this type of regularization did not work even when the contribution of the stabilization term was comparable in magnitude to the contributions from other terms in the TLA codes. We attribute this phenomenon to the specific structure of the unstable modes in the TL equations, which often take the form of strongly anisotropic ridge-like structures (i.e., having a wide spatial spectrum of the SIT component) in the areas of ice convergence. As a consequence, the unstable modes cannot be efficiently damped using isotropic diffusion added to the linearized equations for the σ and/or ice velocity components located in the respective rows of MX.

A straightforward solution is to introduce a spatially inhomogeneous diffusion tensor field (e.g., Yaremchuk and Nechaev, 2013), with local anisotropy derived from the background solution. However, this requires a considerable reduction of δts due to very large diffusion along the ridges. As a simple alternative, we employed Newtonian friction in the TL version of Eqs. (13)–(15), which homogeneously damps the entire spectrum of small perturbations. With this regularization, additional terms εNσi,i=1,2,3 appear inside the square brackets of the linearized Eqs. (13)–(15), where εN is the ratio of δts to the Newtonian damping timescale TN (see Appendix A for details). Numerical experimentation has shown that this approach worked generally well using the Newtonian damping timescale TN of 7δts.

Additional experiments have shown that TN must decrease to 3δts in the case of stronger sea ice convergence in the regions with thick (h>3 m) ice.

Testing the validity of the stabilized TLA codes was done in a way similar to Yaremchuk et al. (2009). The initial conditions for the model thickness and concentration fields C={h(x,0),A(x,0)} were slightly perturbed by the realizations of a random function R (viz. C(x,0)C(x,0)+ϵR(x)C+δC), and the model and its TL version were integrated for t=5δt. After that, the dependence of the normalized difference between the nonlinear solution and its TL approximation was checked by computing the following quantity:

(25) Φ ( ϵ ) = | X ( C + δ C ) - X ( C ) - T L ( δ C ) | | X ( C ) |

, using the initial conditions listed in the argument, and || is the Euclidean norm. As it is evident from Fig. 1, the stabilized version of the TL code is characterized by Φ(ϵ)∝ϵ, while the correct TL code should provide the decay proportional to the square of ϵ. Decay of Φ for the unstable TL code is slightly faster than the stabilized one, but the respective solutions produce very noisy patterns causing much earlier stagnation of the descent process as compared to the stabilized code.

Figure 1TL approximation errors Φ(ε) of the EVP (black), regularized EVP (red) and 1D VP ice model (blue dashed) solutions.


It is noteworthy that the behavior of Φ(ϵ) in similar experiments with the 1D VP sea ice model and corresponding TLA models (Stroh et al., 2019) agreed well with the ϵ2 dependence of the Taylor expansion (blue line in Fig. 1). We speculate that, similar to the 1D case, stability criteria of the 2D VP system may also be weakly affected by the TL transition due to the specific nature of the linearization in the implicit momentum equation solver (see Appendix B). These considerations may indicate an additional attractive feature of the VP rheology in the practical 4D-Var applications. Note also that, to the best of our knowledge, the existing 4D-Var sea ice models are based on VP solvers (MITgcm, Heimbach et al., 2010; NAOSIM, Kauker et al., 2009) and have never reported instability of their TLA codes. Our experiments with a larger number of internal iterations (up to 2000) did not reveal any substantial difference in either nonlinear model solutions or in stability of the TLA models. More detailed descriptions of the TLA codes are given in Appendix A.

In the OSSEs described below, control variables include initial conditions for sea ice velocity, thickness and concentration; the wind stress; and the spatially varying RP fields of P*, e, kT and k2. For other RPs, we utilized constant values adopted from Lemieux et al. (2015, 2016) and listed in Table 1.

2.3 Simulated observations and cost functions

In all OSSEs we used three types of simulated sea ice observations, trying to keep the magnitude of sea ice observational errors close to realistic values.

The first type of data are accurate SIC observations, of which there are currently multiple gridded products based on various remote sensing instruments with different spatial resolutions. After additional preprocessing, these observations are routinely used in data assimilation systems (e.g., GOFS 3.1 DA system of Cummings and Smedstad, 2013), with a nominal spatial resolution of 5 km and regionally low SIC representation errors (5 %, Yaremchuk et al., 2019).

The second data type are SIT observations which contain moderate errors. Currently, the primary source of such data is CryoSat-2, with 1 and 5 km gridded 2 d averaged observations available from the Center for Polar Observations and Modeling (, last access: 1 December 2020). Currently, the error estimates of CryoSat-2 SIT observations range between 0.34 and 0.74 m (Alexandrov et al., 2010; Laxon et al., 2013; Tilling et al., 2018). The recently launched ICESat-2 satellite provides high-resolution (40 m) freeboard estimates (, last access: 1 December 2020) with higher accuracy. So, in the future, novel observational platforms and methods of analysis will likely provide better spatial coverage (i.e., over the entire Arctic) and improved accuracy. In our experiments, we set SIT observation errors to 0.3 m, having in mind future improvements of the SIT data accuracy. A similar error level was adopted by Stroh et al. (2019), which studied RP retrievals in a 1D sea ice model.

Accurate observations of sea ice velocities compose the third data type. An example product is the daily 25 km SIV analysis of various satellite and in situ sources (Tschudi et al., 2019). The respective uncertainties were established at 0.01–0.02 m/s (Schwegmann et al., 2011; Sumata et al., 2015). New methods of sequential SAR image comparison can resolve high-resolution SIV with an accuracy of 0.005 m/s (Komarov and Barber, 2014), suggesting a possibility of high-precision SIV observations in the near future. In the OSSEs reported below, inaccuracy of SIV is set to 0.01–0.025 m/s. Simulated SIC, SIV and SIT observations were derived from the true solution by adding the abovementioned errors with a spatial decorrelation scale of 150 km and a temporal decorrelation scale of 7 d. Taking into account that high-resolution satellite SIC and SIV are currently available on a daily basis over the entire Arctic Ocean and assuming they can be interpolated within each daily time frame, we set the synthetic observations to be available in all the space–time grid points of the model domain.

The standard state-space 4D-Var DA approach of Le Dimet and Talagrand (1986) was utilized: the optimal vector C of control variables was sought to ensure that observations of the model states lie close to assimilated observations within the prescribed time interval (assimilation window), which was set to 3 d in all the experiments. The DA procedure was performed by minimizing the quadratic cost function 𝒥m(C), which included simulated data and regularization (smoothness) terms both characterized by the diagonal error covariance matrices (Appendix A). In addition, we applied bounding constraints on the field values of ice concentration (0A1) and the control fields of the rheological parameters (listed at the bottom of the left column in Table 1).

Table 2List of the performed experiments

Download Print Version | Download XLSX

2.4 OSSEs

The major goal of the conducted OSSEs was to evaluate the feasibility of reconstructing the RPs through assimilation of the SIV, SIC and SIT observations. The first group OSSEs (named KT and K2) analyze the feasibility of optimizing landfast ice parameters k2 and kT. The control vector in these experiments included only parameters k2 and kT, and first-guess initial conditions for SIT and SIC were assumed to be error-free and were not adjusted. Observations of the sea ice velocities, thickness and concentration were assimilated. In the second group of experiments forced by gyre-shaped winds (GYRE-0, GYRE-W), we analyzed the feasibility of optimizing P* and e in the regions with spatially and temporally varying SIT and SIC. In the third group, we analyzed the feasibility of optimizing P* and e in the pack ice zone (PIZ), where SIT varies in space and SIC is equal to 1. We also explore the impact of optimizing P* and e on the hindcast of ice thickness and internal stress distributions. In the second and third groups of the experiments, both the RPs and initial conditions were optimized, and the first-guess sea ice initial conditions were, accordingly, disturbed. A list of OSSEs and their short descriptions are assembled in Table 2. The maximum number of control variables associated with the initial conditions (the number of ice model grid points occupied by the SIT, SIC and SIV fields) was about 9000. As mentioned in Sect. 2.2.1, the RP control fields were defined on coarser (δxp=15δx) grids and bilinearly interpolated on the model grid of the respective OSSEs. Thus, the maximum dimension of the RP control vector never exceeded 36=(75/15+1)(30/15+1) elements, where 75 and 30 are the grid dimensions in Table 1. In all the experiments, we assumed that SIT, SIC and SIV observations were available at all the space–time grid points of the model domain.

3 Optimization of the landfast parameters

3.1 Arching: optimization of kT

Formation of landfast ice in the deep narrow straits and between islands is a well known phenomenon in the Canadian Archipelago and in the Kara Sea (e.g., Lemieux et al., 2016). In the Nares Strait, landfast ice is observed periodically and typically its boundary has an arching shape (e.g., Ryan and Münchow, 2017).

To mimic this phenomenon, the sea ice model was configured in a narrow zonal domain forced by steady zonal wind for 3 d (Fig. 2a). The initial distributions of SIT/SIC were zonally symmetric and SIT/SIC fields were set with values of 2 m∕1.0 and 0.2 m∕0.7 in the western and eastern parts of the domain, respectively, while initial velocities were set to zero (Fig. 2a). Following König Beatty and Holland (2010) and Tremblay and Hakakian (2006), the true value of kT was set to 0.6. Figure 2b shows that after 3 d, the sea ice in the western part of the domain did not drift eastward due to internal tensile strength, which was strong enough to keep sea ice in place, i.e., forming landfast ice in the western part of the domain. In the eastern part where the tensile strength was weaker due to thinner (0.2 m) ice and smaller SIC (0.7), the sea ice moved eastward with typical velocities of about 0.1 m/s forming a polynya between the landfast ice area and thinner sea ice. Due to the impact of the Coriolis force, ice moved slightly southward, forming the polynya along the northern boundary and increasing SIC along the southern boundary.

Figure 2Results of the OSSEs optimizing kT. (a, b) True SIC and SIV with kT=0.6 at t=0 (a) and t=3 d (b). (c, d) The first-guess (c) and optimized (d) SIC and SIV fields at t=3 d. (e, f) Evolution of the normalized cost function for the OSSE with optimization of kT only (solid line) and with the joint optimization of kT, h and A (dashed line) (e). Left part of panel (f) shows optimized kT for the experiment with true first-guess SIC/SIT distributions at t=0. (g, h) Results of the OSSEs optimizing kT with true kT=0.2 and zonal wind of 6 m/s: optimized velocity field and concentration at t=3 d (g); optimized kT (h).


The first-guess solution was forced by the same wind but with kT=0. Figure 2c shows the first-guess state at the end of the assimilation window (t=3 d). It clearly shows that SI moved eastward with a speed of 0.1–0.15 m/s throughout the entire domain, with the sharp boundary between thick and thin ice (Fig. 2a) deteriorating and the 10 km wide polynya developing at the western boundary. Since in this solution ice moves eastward over the entire domain, landfast ice is completely absent due to the absence of tensile strength in the ice (kT=0).

The 4D-Var optimization of only kT (initial distributions of SIC and SIT were not optimized) provides a significant improvement of the SIC and SIT (latter not shown) clearly seen in Fig. 2d. The optimized kT (Fig. 2f) is very close to the true (0.6) value almost everywhere in the western part of the domain, while in the eastern part it is close to zero. Thus, the optimization of kT enabled formation of landfast ice in the western part of the model domain.

Figure 3Results of the k2 optimization with k̃1=0.8. (a, b) True SIC and SIV with k2= 15 at t= 0 and t= 3 d respectively. SIT distribution (meters) is shown by white contours in the left panel. (c, d) The first-guess SIC and SIV with k2= 0 at t= 3 d (c) and (d) optimized SIC and SIV at t= 3 d. (e, f, g) Zonal topography profile (e), evolution of the normalized cost function for the OSSE with optimization k2 (f) and the optimized k2 distribution (g).


Obviously, a similar effect could be achieved with much higher values of kT. To remove this ambiguity, we added an additional term into the cost function which is proportional to the integral of kT2 over the domain. By minimizing the magnitude of kT, we find the minimum value of kT necessary for holding ice in place. Note that optimization of kT was achieved in only four iterations (Fig. 2e), which is a consequence of our sparse-grid representation of RPs in the 4D-Var experiments. If the initial SIT and SIC distributions are not error-free, it is also possible to optimize the initial SIT and SIC in addition to kT. This kind of optimization provides better SIT/SIC hindcast but requires more iterations to find the optimal solution (dashed line in Fig. 2e).

Lemieux et al. (2016) conducted multiple experiments with different values of the kT specified over the entire Arctic Ocean and found that kT should be smaller than 0.6, the value originally proposed by König Beatty and Holland (2010). Because of this, we conducted an additional experiment in which the value of kT was set to 0.2 everywhere and a weaker wind of 6 m/s was specified. Sea ice thickness and concentration were the same.

The optimized solution after 3 d (Fig. 2g) is very similar to the solution in the experiment with kT=0.6, but velocities are smaller and sea ice concentrations in the polynyas are higher due to weaker wind forcing. However, the spatial pattern of the optimized kT distribution is different: it has a clear meridional maxima of 0.18–0.21 between 450 and 600 km and sharply decreases to nearly zero in the other parts of the domain. The largest maximum of the optimized kT values is very close to the true value of kT= 0.2 utilized for this experiment. Note that ice velocity in the entire western part is still equal to zero because the optimized tensile strength is sufficiently strong to keep all ice in place, while wind is not strong enough to deform the sea ice in the western part of the domain. This result suggests that accurate landfast ice modeling can be achieved by specifying nonzero kT only in the key regions, and thus there is no need to specify uniform kT as it was done in the experiments conducted by Lemieux et al. (2016). In operational practice, such arching regions can be identified by analyzing SAR and/or sea surface temperature (SST) images (e.g., Ryan and Münchow, 2017), or from historical sea ice maps available from sea ice data centers.

3.2 Grounding effect: optimization of k2

Grounding on the shallows is another mechanism of landfast ice formation. This kind of landfast ice is typically observed in the Laptev, Chukchi and East Siberian seas and along the northern Alaskan coast (e.g., Lemieux et al., 2015). To mimic this phenomenon, the model was configured in the rectangular 1125×450 km domain (Fig. 3) with zonally varying depth ranging between 3 m at the western boundary and 33 m at the eastern boundary (Fig. 3e). The model was forced by the uniform zonal 10 m/s wind. The true solution was specified as follows. The initial SIC ranged between 0.2 and 1, while the initial SIT was proportional to the SIC and ranged between 0.25 and 2.5 m as shown in Fig. 3a. Initial velocities and the tensile/compressive strength ratio kT were set to zero, while the following true values of the RPs (Lemieux et al., 2016) were used: the critical thickness parameter k̃1= 8, the basal stress parameter α̃b= 20 and the maximum basal stress parameter k2= 15 N/m3. Figure 3b shows that after 3 d the sea ice moved eastward in most of the domain with a typical speed of 0.2–0.4 m/s; only in the southwestern corner was the combination of SIT, SIC and bottom topography sufficient to keep sea ice in place, thus forming a region of grounded landfast ice. Another interesting feature is the elongated polynya visible along the eastern boundary of the landfast ice region approximately 400 km from the western coastline (Fig. 3b). The first-guess solution had the same initial conditions, wind forcing and RP distribution with the exception that k2 was set to zero.

Initial SIT and SIC fields were set to the true values. Despite the perfect initial condition, it is clearly seen that wind moves sea ice eastward with a speed of 0.3–0.4 m/s and forms a polynya along the western boundary (Fig. 3c), which does not exist in the true solution (Fig. 3b). The landfast ice region (i.e., the area with no SIV) in the southwestern corner is completely absent. This polynya separating the landfast ice from the moving ice in the south (bottom of Fig. 3b) is also absent in the first-guess solution.

The variational assimilation of SIV, SIT and SIC observations, targeted at optimization of k2, demonstrated a significant improvement of SIT (not shown) and SIC, clearly seen in Fig. 3d. In particular, the optimized solution includes a landfast ice region and a polynya, which are nearly identical to those in the true solution (Fig. 3b).

The optimized field of k2 is shown in Fig. 3g. The maximum values of k2 are very close to the true k2= 15 N/m3. Note that the true k2 value was specified ad hoc and the grounding effect formally could be reached with smaller values of k2. This is clearly demonstrated in Fig. 3g, where k2 is about 12 N/m3 over the major part of the landfast ice area in true solution (cf. Fig. 3b, g).

Similar to the KT experiment, an additional term penalizing the magnitude of k2 was added to the cost function, and the optimized field of k2 was obtained after a relatively small number of iterations (10–12) (Fig. 3f). Note also that, according to Eqs. (1) and (8), sea ice acceleration is directly proportional to k2, which should invoke faster convergence of the minimization procedure.

To demonstrate the robustness of the k2 optimization with respect to possible variations in k̃1, we conducted a similar experiment with a smaller true k̃1= 2.5. As seen in Fig. 4a and b, the decrease in k̃1 causes a considerable decrease in the landfast area in the southwestern corner. The landfast ice polynya has also moved approximately 100 km closer to the western coastline and is now confined to the shallower region. The first-guess solution for this experiment was the same as in Fig. 3c with k2= 0, i.e., all ice moving eastward. However, after 4D-Var DA and optimization of k2, the reconstructed solution (Fig. 4d) is nearly identical to the true solution (Figure 4b). This optimized k2 (Fig. 4c) has the same spatial structure but a slightly smaller (  14 N/m2) maximum value compared to the experiment where k̃1= 8.

Figure 4Results of k2 optimization with k̃1=2.5. (a, b) True SIC (color) and SIV distributions with k2=15 at t= 0 (a) and t= 3 d (b). SIT distribution (meters) is shown by white contours in the left panel. (c, d) Optimized k2 distribution (c) and optimized SIC and SIV at t= 3 d (d).


The parameterization of grounding landfast ice also includes the critical thickness parameter k̃1, which was kept fixed in the described experiments. According to multiple numerical simulations, the total landfast ice area is more sensitive to variations of k̃1 than k2 (Lemieux et al., 2015) because k̃1 can be interpreted as a scaling coefficient in the definition of the critical ice thickness k̃1hc=Ahb (cf. Eq. 8). It, therefore, acts as a switch which defines the areas of potential landfast ice generation. However, the discontinuity of the Heaviside step function present in Eq. (8) significantly complicates k1 optimization through the gradient-based variational method. Formally, this problem can be regularized (e.g., Nicolsky et al., 2009), but such an approach requires the additional optimization of a regularization parameter, which may ultimately be less computationally efficient in practice. In light of this consideration, we limit our feasibility analysis of landfast ice parameterizations to optimizing kT and k2.

Figure 5True solution in GYRE-0/W OSSEs: (a–b) evolution of the SIT, SIC (white contours) and SIV (black arrows) at t= 0 (left panel) and 3 d; (c–d) evolution of the trace Ptr of the internal stress tensor at t= 0 (left) and 3 d. White arrows show the initial wind stress caused by specified cyclonic wind forcing; (e–f) true distributions of P* (kN/m 2, left panel) and e.


4 Optimization of the ice strength and axes ratio fields

4.1 Cyclonic gyre experiments (GYRE-0/W)

The rheological parameters P* and e are the most important set of parameters responsible for proper sea ice modeling in the deep part of the Arctic Ocean. To evaluate the feasibility of optimizing P* and e, the EVP model was configured for a 2250×900 km rectangular domain with initial true values of the SIC and SIT fields as shown in Fig. 5a. The true values of P* and e varied as shown in Fig. 5e and f within the following ranges: 22.5 kN/m2P*32.5 kN/m2 and 1e 3. These ranges were adopted from various studies (e.g., Hibler and Walsh, 1982; Kreyscher et al., 1997, 2000; Tremblay and Hakakian, 2006; Lemieux et al., 2016). The true wind forcing had a form of a Gaussian-shaped cyclone with a stationary position whose strength gradually increased by 1.5 times during the 3 d assimilation window. The resulting wind stress at t= 0 is shown in Fig. 5c and had a maximum value of 0.7 N/m2. Initial SIV conditions were determined by a 100 min model integration starting from rest, with the all other initial variables and parameters being the same. The initial internal stress was small, but significantly increased after 3 d, under the applied atmospheric forcing. Figure 5c and d show the trace of the stress tensor σIPtr=-trσ/2 at the beginning of the true state integration and after 3 d. The Ptr distribution has a clear maximum near the location with the coordinates (500 km, 500 km), which corresponds to the maximum pressure ( 40 kN/m2) in the sea ice field (e.g., Tremblay and Mysak, 1997). This maximum is due to strong convergence of the relatively thick ( 2.5 m) sea ice in this region. In the eastern part of the domain, Ptr is typically very low due to the divergence of the SIV and considerably thinner ( 0.5–1.5 m) sea ice.

Noisy SIC, SIT and SIV observations were generated by adding spatially and temporally correlated noise (with the correlation scales of 150 km and 7 d) to each of the state variables of the true solution at every time step. The simulated data mimics realistic observations such as those obtained from sources discussed in Sect. 2.3, i.e., they have similar absolute errors to most of the currently available observations. In the experiments, we did not introduce any bias to ice observations since the bias-free observations are a common assumption in existing DA systems.

The magnitudes of the imposed noise correspond to errors of the respective observational data sets, with the amplitudes of 0.05, 0.25–0.35 m, 0.025 m/s and 0.01–0.025 N/m2 for SIC, SIT, SIV and wind stress, respectively. The initial conditions for the first-guess solution were generated in a way similar to the true solution, with slightly larger decorrelation length scales for SIT, SIC and SIV and spatially uniform values of P*= 27.5 kN/m2, e= 2 and true wind forcing.

Figure 6Results of the GYRE-0 OSSE. Suboptimal distribution of SIT, SIC (a) and Ptr0 (c) for t= 3 d after optimization of the initial conditions u0, h0 and A0 using the first-guess values of P*= 27 kN/m2 and e= 2. The deviation norms Su and Str are shown. Panels (b) and (d) are the same as panels (a) and (c) but after additional optimization of P* and e; (e, f) optimized distributions of P* (kN/m2) and e. Black arrows in panels (a) and (b) show the difference between optimized and true SIV.


Despite the exact wind forcing, the first-guess solution differs significantly from the true solution. Similarly to Stroh et al. (2019), the optimization was conducted in three steps. First, we optimized initial SIV, SIT and SIC conditions C0={u0,h0,A0}. Then we sequentially optimized rheological components of the control vector Crh={P*,e} and finally conducted an additional optimization of the full control vector C={C0,Crh}. Note that available SIV, SIT and SIC observations efficiently constrain the respective initial conditions and thus provide a more accurate first guess for the final optimization of the entire control vector. This is important for highly nonlinear optimization problems whose cost functions may have multiple minima.

Figure 6a–d compare the model states after optimization of the initial conditions C0 with fixed P* and e (left panels) and after additional spatial optimization of the P* and e (right panels). The major result of optimizing Crh is an improvement of the SIV fields. The formal quantitative measure of the SIV improvement was evaluated by the function


where Ω is the model domain and index k=1,2 enumerates suboptimal optimization stages: k=1 using the initial conditions control vector C0 only and k=2 employing the full control vector C={C0,Crh}. It was found that Su reduced almost 1.5 times after the additional optimization of the rheological parameters Crh (Fig. 6a, b). Visual comparison of the suboptimal and fully optimized SIC shows a certain improvement after Crh as well. For example, the local minimum of the suboptimal SIC in the region with coordinates 700–1500 and 180 km is about 0.7 (white contours in Fig. 6a), while the fully optimized SIC has a minimum of 0.6 and agrees perfectly with the true SIC distribution (white contours in Fig. 5b).

Figure 7Results of the GYRE-W OSSE: optimized distribution of SIT, SIC (a) and Ptr (b) at t= 3 d. Black/white arrows show the differences between true and optimized SIV (a) and wind stress (b) respectively; (c, d) optimized distributions of P* (kN/m2, left) and e.


In this experiment, we found that optimization of Crh yields only a marginal improvement of the SIT distribution. In particular, SD(hoptk-htrue) decreased from 0.23 m, after optimization of the initial conditions, to 0.2 m, after additional Crh optimization. The minor impact of Crh optimization on the SIT is probably due to relatively high SIT errors and a substantial difference between the first-guess and observed SITs. Another possible reason is that the initial SIT distribution is not well controlled by Crh over the relatively short timescale of the 3 d assimilation window.

As expected, optimization of the rheological parameters Crh provides a major impact on the correction to the internal stress tensor. Figure 6c and d show that standard deviation of the differences between the true and the optimized values of Ptr decreased by  35 %, after the additional RP optimization. Note, also, that the fully optimized Ptr maximum (Fig. 6d) is located in close proximity to the true Ptr maximum (Fig. 5d), while the maximum in the suboptimal Ptr is shifted almost 400 km northward. Comparing optimized and true P* fields demonstrates agreement almost everywhere, with the only exception observed in the southwest corner of the domain (Fig. 6e). This is caused by the substantial SIV divergence (Fig. 6a), which diminishes the role of rheological forcing in the region.

The reconstruction of e is not as accurate as P*; the optimized e ranges between 1.2 and 2.8, while the respective true values are between 1 and 3. However, the spatial locations of the extrema are in strong agreement. The exception is the local minimum in the southwestern corner, where both e and P* disagree with true values and are close to their first-guess values. Note that significant improvement of Ptr and the internal stress components discussed above are directly related to the more accurate optimization of the RP and SIV fields because both P* and e control the structure of the stress tensor in Eq. (2).

To analyze the impact of wind forcing inaccuracy, we conducted an additional experiment where the center of the cyclonic disturbance was displaced 90 km westward, mimicking a systematic error in the hypothetical atmospheric forecast. The results obtained after full optimization of the control vector C={C0,Crh} are shown in Fig. 7. It is worthy to note that the inaccurate position of the cyclone causes significant errors (up to 0.2 N/m2, or  25%) in the wind stress forcing in the central part of the domain (Fig. 7b). As a result, the optimized SIV fields have essential ( 0.1 m/s) errors (Fig. 7a) and the integral measure of the SIV inaccuracy Su increased five times up to 0.64 m/s.

At the same time, degradation of the SIT retrieval was not as significant, with SD(hoptk-htrue) increasing up to 0.25 m, i.e., by only 25 % as compared to the previous experiment with exact wind forcing. Similarly, the optimized SIC distribution remained largely unchanged. The integral quality of the reconstruction of Ptr is 3.5 kN, i.e., about 40 %–50 % worse than in the experiment with exact forcing, but it is important to note that maximum of Ptr is remains in very good agreement with the true solution.

Although inaccurate wind forcing has a profound impact on the accuracy of P* and e retrievals, there is still an essential level of similarity between the reconstructed and true rheological fields. For example, spatial distribution of the optimized P* still has its maxima in the western and eastern parts of the region and a minimum in the center of the domain (Fig. 7c), while the minimum of the e, in the western part, demonstrates a certain agreement with true e distribution (Fig. 5e, f). Note that inaccurate wind forcing affects the accuracy of P* retrievals to a lesser degree than that of e.


In both experiments, described in the previous section, the initial true SIC was rather close to 1 and decreased below 0.8 in some regions after 3 d of integration (Fig. 5a, b). Due to the exponential dependence of P on 1−A (Eq. 7), the internal stress decreases exp (4)≈50 times and therefore has a minor rheological impact on the sea ice dynamics in these regions. At the same time in winter, most of the Arctic Ocean is almost completely covered by pack sea ice, with SIC ranging between 0.98 and 1. To mimic these conditions, we conducted another OSSE with spatially and temporally invariant sea ice concentration A=1. Numerically, this was achieved by removing the advection equation from Eq. (4) and removing initial A0 from the control vector C0. Note that setting A=1 can be interpreted as the introduction of a very fast cooling, which immediately removes areas with A<1.

The model domain, initial SIT, and P* and e distributions were the same as in GYRE-0/W OSSEs. However, unlike the cyclonic wind from the previous experiment, the applied atmospheric forcing is a 16 m/s eastward wind at the western boundary which reverses in zonal direction across the breadth of the domain (Fig. 8c). In time, the wind speed was linearly increased up to 20 m/s by the end of the DA window. The resulting wind stress at t= 3 d (Fig. 7c) has a maximum amplitude of about 0.5 N/m2. This relatively strong wind corresponds to the category of strong Arctic cyclones, a rather typical phenomenon in the Arctic Ocean, which may persist for periods of up to 2 weeks (Simmonds and Rudeva, 2012).

Figure 8True solution in PIZ OSSEs: (a–b) Evolution of the SIT and SIV (black arrows) at t= 0 (left) and 3 d; (c–d) evolution of the Ptr of the internal stress tensor. White arrows show the wind stress forcing at t=0. The distributions of P* (kN/m2) and e are the same as for GYRE-0/W experiments (Fig. 5e, f).


The temporal evolution of the true SIT and SIV is shown in Fig. 8a and b. Under the relatively strong applied wind forcing, sea ice converges and SIT increases almost everywhere, with the exception of the narrow bands along the western and northern boundaries, caused by the joint effect of the Coriolis force and coastal boundary conditions. The true SIV has a maximum of about 0.4 m/s. The distribution of Ptr, shown in the lower panels of Fig. 8, has two clear maxima in the south, with the magnitudes of about 70 and 50 kN/m2. Note that due to ice convergence causing SIT growth both Ptr maxima increase in magnitude and slightly ( 60 km) move towards each other.

The first-guess initial SIT/SIV conditions and data for the PIZ OSSE were derived from the true solution in a similar way as for the GYRE-0 OSSE. Similarly, we specify the first guess with P*= 27.5 kN/m2, e= 2 and the exact wind forcing. The optimization was conducted in three steps: by first optimizing C0, then optimizing Crh, and, finally, simultaneously optimizing C0 and Crh.

Figure 9a and c show SIT and the difference between optimized and true SIV at t= 3 d after optimization of the initial conditions only. Interestingly, despite less spatial variations in wind forcing, optimization of the initial conditions does not allow for accurate reconstruction of SIV as in the GYRE-0/W OSSEs. The maximum errors in the eastern part of the domain are about 0.1 m/s, being comparable in magnitude with regional velocities. The relative accuracy of the SIV reconstruction, in the western part, is slightly better but is still considerably worse than in the previous OSSEs.

Figure 9Results of the PIZ OSSE: suboptimal distributions of SIT (a) and Ptr (c) at t= 3 d after optimizing the initial conditions u0,h0,A0 and the first-guess values of P*= 27.5 kN/m2 and e= 2. The discrepancies of Su and Str with the true solutions are shown. Black arrows show the difference between optimized and true SIV for each optimization stage. Panels (b) and (d) are the same as panels (a) and (c) but after additional optimization of P* and e. Panels (e) and (f) show optimized distribution of P* (kN/m2) and e. Panels (g) and (f) show optimized distribution of P* (kN/m2) and e for the experiment with moderate wind of 10 m/s.


The suboptimal distribution of Ptr (Fig. 9c) differs significantly from the true Ptr (Fig. 8d), both quantitatively and qualitatively. For example, SD(Ptr1(opt)-Ptr(true)) is 11.6 kN/m2, or about 30 %–35 % of the domain-average value of Ptr. The qualitative difference is probably more important because the suboptimal Ptr distribution fails to provide the two maxima discussed above. Instead, the suboptimal distribution of Ptr has only one maximum located in the center (Fig. 8d).

Additional optimization of the rheological parameters, Crh, significantly improves the reconstructed SIV practically everywhere, with the formal measure of the uncertainty Su decreasing by almost one-half from 56 to 30 m/s (Fig. 9b). Similar improvements are visible in the Ptr distribution (Fig. 9d). The SD[Ptr(opt)−Ptr(true)] decreased to 7.4 kN/m2 (by 40 %), and the fully optimized Ptr has two maxima as in the true Ptr distribution (Fig. 8d).

The optimized P* and e are shown in Fig. 9e and f. In the eastern part of the model domain, the reconstructed P* and e almost perfectly agree with the true distributions of P* and e, while there is some difference between optimized and true rheology in the western part. This is probably due to offshore sea ice transport which creates ice divergence along the rigid western boundary. As a result, the impact of rheology on ice dynamics becomes less significant here, and rheological parameters are harder to recover. There is also some quantitative difference between optimized and true e in the central part of the model domain, but, qualitatively, the reconstructed field of e has all the features of the true distribution.

As mentioned above, PIZ experiments were conducted under a relatively high wind forcing. To assess the possibility of reconstructing P* and e under weaker winds, we conducted an additional experiment with a typical wind of about 10 m/s and a maximum wind stress about 0.15 N/m2. The results, shown in Fig. 9g and h, are similar to those of the experiments using stronger wind. In particular, both P* and e are reconstructed better in the eastern part of the domain and somewhat worse in the western part due to ice divergence near the western boundary. Compared to the stronger wind case, the quality of P* reconstruction is diminished under weaker winds because lower wind stress generates smaller ice velocities, amplifying the impact of observation errors.

Note, however, that most of the inaccuracies in the reconstructed P* distribution are observed in the region with relatively thick (2.0–2.6 m) ice located meridionally between 700 and 1200 km in the meridional direction. Under moderately weak wind, therefore, it is still possible to accurately reconstruct RPs in thin ice regions. Obviously, RPs are not reconstructible in regions where winds are too weak to influence ice motion. In this case, optimized P* may freely take any value above a critical threshold determined by ice strength.

5 Summary and discussion

The presented study continues our previous efforts (Stroh et al., 2019) and addresses the feasibility of retrieving spatially varying RPs through the 4D-Var assimilation of the satellite observations of SIV, SIC and SIT in two dimensions. To perform this analysis, we developed TLA codes with respect to all rheological parameters (except k1), initial conditions and wind forcing for a single-category sea ice model recently proposed by Lemieux et al. (2016). The dynamical core of this model is based on the conventional formulation of the EVP rheology (Hunke and Dukowicz, 1997; Hunke and Lipscomb, 2008) and parameterizations of the grounding and arching landfast ice recently proposed by Lemieux et al. (2015, 2016) and König Beatty and Holland (2010). The model was configured in multiple rectangular domains and included several simplifications. In particular, we constrained ourselves to a single ice category, utilized a relatively simple but still widely used parameterization for the internal pressure (Eq. 7), and employed a simplified nonlinear Lax–Wendroff scheme for ice advection. We adopted these simplifications to reduce the complexity of the TLA codes, and this had negligible impact on the results at the 3 d timescale of conducted experiments.

It was found that TLA models for the EVP solver are unstable for the regions with high (>0.9) sea ice concentration and require stabilization. The standard stabilization technique, through the additional diffusion (Hoteit et al., 2005), widely used in the OGCM inverse modeling, was found to be inefficient, but a simpler stabilization, based on Newtonian friction, appeared to work well. Analysis of the TL approximation accuracy has shown that Newtonian stabilization has errors similar to the ones observed in the case of diffusion-based stabilization, and thus the Newtonian scheme can be successfully used in sea ice models based on the EVP solvers. In the last decade, several modifications of the EVP solvers were proposed to improve the convergence (e.g., Bouillon et al., 2013; Kimmritz et al., 2016; Koldunov et al., 2019). We assume that these ideas could benefit the development of the respective TLA models. In particular, the Newtonian damping coefficient could be adjusted locally to account for the background sea ice pressure/tensile strength and thus provide a more efficient stabilization of the EVP TLA codes.

On the other hand, simple analysis (Fig. 1, dashed line) indicates that the TLA model for the 1D VP sea ice solver (Stroh et al., 2019) is stable, provided that the forward model satisfies the stability criterion (which is always the case). Similar indications were obtained by Heimbach et al. (2010) and NAOSIM (Kauker et al., 2009), which did not observe instabilities in the TLA codes of the sea ice models with VP rheologies. We may, therefore, conclude that numerical models with VP rheology (e.g., Hibler, 1979; Lemieux et al., 2008) do not require additional stabilization of the TLA codes and are formally more suitable for development of sea ice 4D-Var DA algorithms. Note, also, that due to their stability, VP TLA codes can be easily incorporated into the parent sea ice models to provide improved Jacobian-free Krylov solvers for VP rheology. This approach was employed in the ocean model by Nechaev et al. (2005), where the internal matrix-free bi-conjugate gradient (BiCG) solver was constructed using the model's TLA codes in the implicit–explicit (IMEX) algorithm, where the preconditioned flexible general minimum residual (GMRES) algorithm was implemented in the Jacobian-free mode (Lemieux et al., 2014; Losch et al., 2014).

In a comprehensive series of OSSEs with a simplified EVP sea ice model, it was demonstrated that Newtonian stabilization of the TLA codes allows for a reasonable reconstruction of the RPs. The numerical experiments included two groups of 4D-Var DA experiments.

First, we analyzed the possibility of optimizing the RPs in two different landfast ice parameterization schemes incorporated in the CICE model (Hunke et al., 2010). In the current CICE6 version, all landfast ice parameters are treated as spatially uniform variables, which degrades the accuracy of landfast ice simulations in various parts of the Arctic Ocean (Lemieux et al., 2015), especially in the shallow seas and narrow straits of the Canadian Archipelago. In this study, these parameters were specified by spatially variable functions with a reduced number (10–36, Sect. 2.4) of free parameters that were optimized to fit surface observations. The conducted OSSEs demonstrate that spatially varying landfast ice parameters k2 and kT, which are responsible for grounding and arching phenomena, can be optimized at a relatively low computational cost (5–12 iterations on a sparse grid). We found that the impact of spatially uniform kT could be achieved by specifying kT only along a relatively narrow landfast ice boundary (Fig. 2h), which works as a barrier to prevent ice drift under moderate wind conditions. This observation suggests that parameterizing landfast ice by spatially uniform kT can be reduced to specifying nonzero kT only within these localized regions of the domain. Interestingly, the optimization of both k2 and k1 requires a very small ( 10) number of iterations, which suggests the possibility of efficiently incorporating their optimization into a pan-Arctic operational model only in regions of potential landfast ice arching or grounding phenomenon.

Taking into account that landfast ice typically forms in shallow/coastal areas and in narrow straits, the landfast phenomena can be controlled more efficiently by adding to the control the critical thickness parameter k1, which confines k2 variability to shallow areas and narrow straits areas in a high-resolution pan-Arctic sea ice model. It should be noted that retrieval of k1 is not straightforward because of nondifferentiability of the grounding parameterization scheme (Eq. 8). Optimizing k1 in sea ice models requires further development including more robust constrained minimization tools, such as genetic (Goldberg, 1989) or very fast simulation annealing (Ingber, 1989) algorithms. Another approach is to use parametric regularization of k1, similar to the one utilized by Nicolsky et al. (2009). In this case, TLA models require additional computation of the derivative with respect to the regularization parameter specifying a smooth approximation (e.g, Lemieux and Tremblay, 2009) of the Heaviside function in Eq. (8). We are currently implementing this algorithm for the 2D VP solver.

In the second group of OSSEs, we analyzed the possibility of reconstructing spatially varying sea ice strength P* and ellipse axes ratio e distributions. The respective OSSEs employed simulated observations of SIV and SIC characterized by the root-mean-square errors of 0.025 m/s and 0.05, respectively, while SIT observations were simulated with uncertainties of 0.3 m, which is about 2 times smaller than the errors of the CryoSat-2 observations. The OSSEs also assessed sensitivity of the results with respect to systematic errors in wind forcing.

The OSSEs with accurate (exact) wind forcing demonstrated the feasibility of relatively accurate reconstruction of the P* distribution and less successful, but still reasonable, reconstruction of the axes ratio distribution. Similar results were recently obtained by Stroh et al. (2019), who used a 1D sea ice model featuring VP parameterization to find a higher impact of spatially varying P* than e on the DA quality. We also observed that regions of less accurate P* reconstructions are typically colocated in the regions of strong sea ice divergence and/or SIC concentrations below 0.8, i.e., with the regions where rheology plays a lesser role in sea ice dynamics.

We also found that additional optimization of P* and e (after optimizing the initial state of sea ice) provides a slightly more accurate reconstruction of the SIT and SIC distributions and a significant improvement of the SIV and Ptr fields. Accurate forecasting of Ptr is especially important for maritime use, as it allows vessels to avoid regions with excessive compressive stress.

The OSSE with strong wind convergence in the pack ice (A= 1) demonstrated even better quality reconstructions of e, and especially P*. This can be attributed to the stronger role of internal stress on sea ice dynamics in pack ice. Similarly to the OSSEs with cyclonic wind pattern, we found that additional optimization of RPs provides smaller improvements in the SIT and SIC hindcasts as compared to SIV and Ptr. OSSEs with weaker winds ( 0.0–0.15 N/m2) demonstrated a slight reduction in the accuracy of reconstructed P* and e. This is because weaker winds generate smaller changes in the sea ice state, and observation errors contribute more to the results of assimilation. Note that in the limiting case of zero (or very weak) wind and thick ice, the optimized P* is unconstrained and may take any value sufficient to keep ice in place.

Our 4D-Var applications utilized a relatively small assimilation window (3 d) with an eye towards improving short-term sea ice forecasts in the ice pack and ice edge zones. In these regions, ice rheology typically changes at the timescales of several days (e.g., Panteleev et al., 2019) due to variations in the dynamic and/or thermodynamic forcing. In particular, wind variability may cause profound changes in the ice edge. Meanwhile in the ice pack zones, wind forcing is the major driving factor of polynya formation and ridging processes, where rheological forces become important.

Experiments with cyclonic and zonal wind emphasized the importance of having accurate prior estimates of wind forcing: since sea ice dynamics is significantly controlled by winds (e.g., Thorndike and Colony, 1982), it is hard to expect a reasonable quality 4D-Var reconstruction derived from a model driven by incorrect winds. In this case, both the model simulation and the 4D-Var results will be inaccurate. However, a properly formulated 4D-Var approach may still adjust wind forcing through the assimilation of the SIV/SIC/SIT observations (e.g., Stroh et al., 2019).

Finally, sea ice observations are typically available on a regular (daily) basis and are expected to gain spatial coverage and accuracy in the near future for the SIC and SIV components that are directly observable from space, while SIT observations are still lacking accuracy due to the complex structure and uncertainties in its observation process. Ongoing developments in data acquisition and preprocessing will result in improved sea ice state observability, and this work has demonstrated the feasibility of reconstructing spatially varying RP fields on the basis of these data.

In the present study we utilized realistic observational errors for SIV and SIC, while SIT errors were somewhat smaller than the accuracy of currently available satellite observations. An additional set of experiments with more realistic SIT errors reveals their stronger impact on the reconstruction quality of P* and e, while the reconstruction accuracy of the landfast ice parameters k2 and kT remained virtually unchanged. The incoming satellite platforms (e.g., ICESat-2,, last access: 1 December 2020) with a better SIT observation capability may deliver sufficiently dense and accurate SIT observations required for reasonably accurate estimation of P* and e in the internal regions of the Arctic Ocean.

Analysis of the potential impact of new observations, as well as more realistic inversions employing more complex rheological hypotheses (e.g., the Maxwell elasto-brittle rheology of Dansereau et al., 2016), may be within the focus of our studies in the near future.

Appendix A: Tangent linear numerics

The TL code was obtained by the subtracting a solution of the numerical model X from the evolution equations of the perturbed state, X+δX, and keeping only linear terms in the expansions of all the nonlinearities. Note that the easiest way to conduct this formal procedure is to apply a tangent liner and adjoint model compiler. In particular, this approach was used for the development of the MITgcm and NAOSIM sea ice models (e.g., Kauker et al., 2009) and basal sea ice model (Goldberg and Heimbach, 2013). Because of the availability of multiple automatic tools (e.g., TAMC, Giering and Kaminski, 1998,, last access: 1 December 2020; OpenAD,, last access: 1 December 2020), we briefly outline the major details of the development of the EVP TLA models.

Perturbations of the auxiliary functions ϵ˙1,2,3,m,ζ,Pp, P and Cb of X are given by


Hereinafter, all the terms with variations of the RPs and atmospheric forcing and variations of the quantities that may contain those are underlined. Taking Eqs. (A1)–(A7) into account, the TL equations of Eqs. (13)–(19) are given by


The gradients with respect to RPs and atmospheric forcing control variables are given by


It is important to note that stability of the adjoint model, being strictly related to the stability of the TL model, contains the Newtonian damping given by the boldfaced terms εNδσ1,2,3 in Eqs. (A8)–(A10). The adjoint model is integrated backward in time and requires the solution of the forward model for each time step of the backward subcycling procedure. This can be achieved through either recalculating the forward solution or storing it, which includes storing all the intermediate states of the subcycling procedure, as well as space–time coordinates of the switches in the Heaviside functions. We elected the second option. Note that, formally, implementation of the 4D-Var DA approach does not require the TL model (see Eqs. 2224). However, the development of the TL model cannot be avoided because of the necessity to check its tangent linear property and validity of the respective Lagrangian identities for debugging the transposition procedure of the TL system matrix MX whose action on the state perturbations is represented by Eqs. (A8)–(A14).

Appendix B: On the stability of TL models with VP rheology

The VP system of equations is obtained by eliminating the time derivative in Eq. (2), resolving the remainder with respect to σ,

(B1) σ = 1 e 2 P - I 1 - e 2 2 tr P ,

and substituting (B1) into the momentum equation:


where we used the notation P for the tensor in the rhs of Eq. (B1):

(B5) P = P p ( 1 + k T ) ϵ ˙ - Δ 2 ( 1 - k T ) I / Δ .

The nonlinear system of VP Eqs. (B1)–(B4) is solved in two stages. First, Eq. (B2) is stepped forward using a semi-implicit scheme, which takes the fields of P and Δ from the previous time step, treating only ϵ˙ semi-implicitly. On the second stage, Eqs. (B3)–(B4) are explicitly stepped forward in time using the velocity fields from the first stage (Hibler, 1979; Lemieux et al., 2008). In application to TLA modeling, such a procedure would not provide the exact derivative of the nonlinear scheme with respect to u because Eq. (B2) is not linearized with respect to the variations of Δ(u) and will therefore exhibit behavior similar to the behavior of the regularized code (Fig. 1). Moreover, the exact TLA code of the VP rheology is intrinsically unstable in the regions of ice divergence (e.g., Gray and Killworth, 1995), especially in the areas where tensile stresses associated with arching effects are important.

However, one may expect reasonable performance of the abovementioned incomplete linearization of the VP model in the 4D-Var applications as soon as the stability criterion of the nonlinear model (B1)–(B4) is satisfied. We performed an additional experiment with 1D VP forward and TL models using the modified procedure featuring 10 applications of the GMRES solver and 10 Δ updates on every time step (Lemieux et al., 2008) and did not observe any instabilities in the TL model. In this respect, we may conclude that VP rheology is less susceptible to the TLA instabilities than EVP, which requires introduction of the additional stabilization terms in the TLA code of the stress tensor evolution equation.

Data availability

All inputs to the model experiments are from idealized conditions as specified in the paper. The numerical scheme is provided in large detail in the Appendix and can be easily replicated. Model results in this study are securely stored at the Navy DSRC archive server and can be accessed after obtaining an account at the facility. The corresponding author can be contacted for information to access the archived data once an account has been established.

Author contributions

All authors provided substantial contribution to the models' development, interpretation of the results and writing the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


This research was funded in part under the Naval Research Laboratory’s 6.2 project Modeling Arctic Landfast Ice (program element 62435N) and the 6.1 project Optimization Rheology and Advancing Sea Ice Forecast (program element 61153N). Oceana Francis was supported by the Coastal Hydraulics Engineering Resilience (CHER) Lab, the Civil and Environmental Engineering Department, and the Sea Grant College Program at the University of Hawai`i at Mānoa. Special thanks to the three anonymous reviewers for their comments, which have significantly improved this paper.

Financial support

This research has been supported by the Naval Research Laboratory’s 6.2 project Modeling Arctic Landfast Ice (program element 62435N) and the 6.1 project Optimization Rheology and Advancing Sea Ice Forecast (program element 61153N) funded by the Office of Naval Research.

Review statement

This paper was edited by Christian Haas and reviewed by three anonymous referees.


Alexandrov, V., Sandven, S., Wahlin, J., and Johannessen, O. M.: The relation between sea ice thickness and freeboard in the Arctic, The Cryosphere, 4, 373–380,, 2010. 

Anderson, D. L. and Weeks, W. F.: A theoretical analysis of sea ice-strength, EOS T. Am. Geophys. Un., 30, 632–640, 1958. 

Arnold Jr., C. P. and Dey, C. H.: Observing System Simulation experiments: Past, Present and Future, B. Am. Meteorol. Soc., 67, 687–695,, 1986. 

Bouillon, S., Fichefet, T., Legat, V., and Madec, G.: The elastic–viscous–plastic method revisited, Ocean Model., 71, 2–12, 2013. 

Cummings, J. A. and Smedstad, O. M.: Variational data assimilation for the global ocean, in: Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications, Springer-Verlag, Berlin Heidelberg, 303–343,, 2013. 

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elasto-brittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359,, 2016. 

Errico, R. M.: What is an Adjoint model?, B. Am. Meteorol. Soc., 78, 2577–2591, 1997. 

Giering, R. and Kaminski, T.: Recipes for adjoint code construction, ACM T. Math. Software, 24, 437–474, 1998. 

Gilbert, J. C. and Lemaréchal, C.: Some numerical experiments with variable-storage quasi-Newton algorithms, Math. Progr., 45, 407–435,, 1989. 

Goldberg, D. E.: Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading, MA, 432 pp., 1989. 

Goldberg, D. N. and Heimbach, P.: Parameter and state estimation with a time-dependent adjoint marine ice sheet model, The Cryosphere, 7, 1659–1678,, 2013. 

Gray, J. M. N. T. and Killworth, P. D.: Stability of viscous-plastic sea ice rheology, J. Phys. Oceanogr., 25, 971–978, 1995. 

Harder, M. and Fischer, H.: Sea ice dynamics in the weddell sea simulated with an optimized model, J. Geophys. Res.-Oceans, 104, 11151–11162, 1999. 

Heimbach, P., Menemenlis, D., Losch, M., Campin, J. M., and Hill, C.: On the formulation of sea-ice models. Part 2: Lessons from multi-year adjoint sea ice export sensitivities through the Canadian Arctic Archipelago, Ocean Model., 33, 145–158,, 2010. 

Hibler, W.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. 

Hibler, W. D. and Walsh, J. E.: On modeling seasonal and interannual fluctuations of arctic sea ice, J. Phys. Oceanogr., 12, 1514–1523, 1982. 

Hoteit, I., Cornuelle, B., Kohl, A., and Stammer, D.: Treating strong adjoint sensitivities in tropical eddy-permitting variational data assimilation, Q. J. Roy. Meteor. Soc., 131, 3659–3682, 2005. 

Hunke, E. and Dukowicz, J.: An elastic-viscous-plasticmodel for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997. 

Hunke, E. C.: Viscous-plastic sea ice dynamics with the EVP model: Linearization issues, J. Comput. Phys., 170, 18–38, 2001. 

Hunke, E. C. and Lipscomb, W. H.: CICE: The Los Alamos sea ice model. Documentation and software user's manual version 4.0, Tech. Rep. LA-CC-06-012, Los Alamos Natl. Lab., Los Alamos, NM, 2008. 

Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos Sea Ice Model documentation and software users manual version 4.1, Tech. Rep. LA-CC-06-012, T-3 Fluid Dynamics Group, Los Alamos Natl. Lab., Los Alamos, NM, 2010. 

Ingber, L.: Very fast simulated re-annealing, Math. Comput. Model., 12, 967–973, 1989. 

Juricke, S., Lemke, P., Timmermann, R., and Rackow, T.: Effects of stochastic ice strength perturbation on Arctic finite element sea ice modeling, J. Climate, 26, 3785–3802,, 2013. 

Kauker, F., Kaminski, T., Karcher, M., Giering, R. , Gerdes, R. and Voßbeck, M.: Adjoint analysis of the 2007 all time Arctic sea-ice minimum, Geophys. Res. Lett., 36, L03707,, 2009. 

Kimmritz, M., Danilov, S., and Losch, M.: The adaptive EVP method for solving the sea ice momentum equation, Ocean Model., 101, 59–67,, 2016. 

Koldunov, N. V., Aizinger, V., Rakowsky, N., Scholz, P., Sidorenko, D., Danilov, S., and Jung, T.: Scalability and some optimization of the Finite-volumE Sea ice–Ocean Model, Version 2.0 (FESOM2), Geosci. Model Dev., 12, 3991–4012,, 2019. 

Komarov, A. S. and Barber, D. G.: Sea ice motion tracking from sequential dual-polarization RADARSAT-2 images, IEEE T. Geos. Remote S., 52, 121–136, 2014. 

König Beatty, C. and Holland, D. M.: Modeling landfast sea ice by adding tensile strength, J. Phys. Oceanogr., 40, 185–198,, 2010. 

Kreyscher, M., Harder, M., and Lemke, P.: First results of the Sea-Ice Model Intercomparison Project (SIMIP), Ann. Glaciol., 25, 8–11, 1997. 

Kreyscher, M., Harder, M., Lemke,P., and Flato, G. M.: Results of the sea ice model intercomparison project: Evaluation of sea ice rheology schemes for use in climate simulations, J. Geophys. Res.-Oceans, 105, 11299–11320, 2000. 

Laxon, S. W., Giles, K. A., Ridout, A. L., Wingham, D. J., Willatt, R., Cullen, R., Kwok, R., Schweiger, A., Zhang, J., Haas, C., Hendricks, S., Krishfield, R., Kurtz, N., Farrell, S., and Davidson, M.: CryoSat‐2 estimates of Arctic sea ice thickness and volume Geophys. Res Lett., 40, 732–737,, 2013. 

Le Dimet, F. X.: A general formalism of variational analysis, CIMMS report 22, Cooperative Institute for Mesoscale Meteorological Studies, Norman, OK, USA, 34 pp., 1982. 

Le Dimet, F. X. and Talagrand, O.: Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects, Tellus, 38A, 97–110, 1986. 

Lemieux, J.-F. and Tremblay, B.: Numerical convergence of viscous‐plastic sea ice models, J. Geophys. Res.-Oceans, 114, C05009,, 2009. 

Lemieux, J.-F., Tremblay, B., Thomas, S., Sedlacek, J., and Mysak, L. A.: Using the preconditioned generalized minimum residual (GMRES) method to solve the sea-ice momentum equation, J. Geophys. Res.-Oceans, 113, C10004,, 2008. 

Lemieux, J.-F., Knoll, D., Tremblay, B., Holland, D., and Losch, M.: A comparison of the Jacobian-free Newton-Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: a serial algorithm study, J. Comp. Phys., 231, 5926–5944, 2012. 

Lemieux, J.-F., Knoll, D., Losch, M., and Girard, C.: A second-order accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics, J. Comput. Phys., 263, 375–392, 2014. 

Lemieux, J.-F., Tremblay, L. B., Dupont, F., Plante, M., Smith, G. C., and Dumont, D.: A basal stress parameterization for modeling landfast ice, J. Geophys. Res.-Oceans, 120, 3157–3173,, 2015. 

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.-Oceans, 121 , 7354–7368, 2016. 

Lewis, J. M. and Derber, J. C.: The use of adjoint equations to solve a variational adjustment problem with advective constraints, Tellus A, 37, 309–322,, 1985. 

Losch, M., Fuchs, A., Lemieux, J. F., and Vanselow, A.: A parallel Jacobian-free Newton-Krylov solver for a coupled sea ice ocean model, J. Comp. Phys., 257, 901–911, 2014. 

Massonnet, F., Goosse, H., Fichefet, T., and Counillon, F.: Calibration of sea ice dynamic parameters in an ocean-sea ice model using an ensemble Kalman filter, J. Geophys. Res.-Oceans, 119, 4168–4184,, 2014. 

Massonnet, F., Fichefet, T., and Goosse, H.: Prospects for improved seasonal Arctic sea ice predictions from multivariate data assimilation, Ocean Model., 88, 16–25, 2015. 

Metzger, E. J., Smedstad, O. M., Thoppil, P. G., Hurlburt, H. E., Cummings, J. A., Wallcraft, A. J., Zamudio, L., Franklin, D. S., Posey, P. G., Phelps, M. W., Hogan, P. J., Bub, F. L., and DeHaan, C. J.: US Navy operational global ocean and Arctic ice prediction systems, Oceanography, 27,32–43,, 2014. 

Miller, P. A., Laxon, S. W., Feltham, D. L., and Cresswell, D. J.: Optimization of a sea ice model using basinwide observations of Arctic sea ice thickness, extent, and velocity, J. Climate, 19, 1090–1108, 2006. 

Nechaev, D. A., Panteleev, G., and Yaremchuk, M.: Reconstruction of the Circulation In Limited Regions of an Ocean With Open Boundaries: Climatic Circulation In the Tsushima Strait, Oceanology, 45, 761–780, 2005. 

Nguyen, A. T., Menemenlis, D., and Kwok, R.: Arctic ice-ocean simulation with optimized model parameters: Approach and assessment, J. Geophys. Res., 116, C04025,, 2011. 

Nichols, N. K.: Data Assimilation: Aims and Basic Concepts, in: Data Assimilation for the Earth System, NATO Science Series (Series IV: Earth and Environmental Sciences), edited by: Swinbank, R., Shutyaev, V., and Lahoz, W. A., Springer, Dordrecht, 2003. 

Nichols, N. K.: Mathematical concepts of data assimilation, in: Data assimilation: making sense of observations, edited by: Lahoz, W., Khattatov, B., and Menard, R., Springer, Dordrecht, 2010. 

Nicolsky, D. J, Romanovsky, V. E., and Panteleev, G. G.: Estimation of soil thermal properties using in-situ temperature measurements in the active layer and permafrost, Cold Reg. Sci. Technol., 55 120–129, 2009. 

Nitta, T.: Some analyses of observing systems simulation experiments in relation to the First GARP Global Experiment, GARP Working Group on Numerical Experimentation, Report No 10, 1–35, Plan for U.S. Participation in the Global Atmospheric Research Program, National Academy of Sciences, Washington, DC, USA, 1969. 

Panteleev, G., Rogers, W. E., Yaremchuk, M., Shen, H., Rainville, L., and Grout, J.: Floe Size Mapping from Satellite SAR Images and Icewatch Observations in the Beaufort Sea during Autumn 2015, Tech. Rep. NRL/MR/7322-19-9903, Naval Researh Laboratory, Stennis Space Center, Mississippi, USA, 2019. 

Penenko, V. V.: Methods of Numerical Simulation of Atmospheric processes, Gidrometeoizdat, Lenigrad, 350 pp., 1981. 

Posey, P. G., Metzger, E. J., Wallcraft, A. J., Preller, R. H., Smedstad, O. M., and Phelps, M. W.: Validation of the 1/128 Arctic Cap Nowcast/Forecast System (ACNFS), Tech. Rep. NRL/MR/7320-10-9287, Naval Res. Lab., Stennis Space Center, Mississippi, USA, 2010. 

Ryan, P. A. and Münchow, A.: Sea ice draft observations in Nares Strait from 2003 to 2012, J. Geophys. Res., 122, 3057–3080,, 2017. 

Schwegmann, S., Haas, C., Fowler, C., and Gerdes, R.: A comparison of satellite-derived sea-ice motion with drifting-buoy data in the Weddell Sea, Antarctica, Ann. Glaciol., 52, 103–110, 2011. 

Simmonds, I. and Rudeva, I.: The great Arctic cyclone of August 2012, Geophys. Res. Lett., 39, L23709,, 2012. 

Stroh, J. N, Panteleev, G., Yaremchuk, M., Francis, O., and Allard, R.: Toward optimization of rheology in sea ice models through data assimilation, J. Atm. Oceanic Tech., 36, 2365–2382,, 2019. 

Sumata, H., Kwok, R., Gerdes, R., Kauker, F., and Karcher, M.: Uncertainty of arctic summer ice drift assessed by high-resolution SAR data, J. Geophys. Res.-Oceans, 120, 5285–5301, 2015. 

Sumata, H., Kauker, F., Karcher, M., and Gerdes, R.: Simultaneous Parameter Optimization of an Arctic Sea Ice–Ocean Model by a Genetic Algorithm, Mon. Weather Rev., 147, 1899–1926,, 2019. 

Thorndike, A. S. and Colony, R.: Sea ice motion in response to geostrophic winds, J. Geophys. Res., 87, 5845–5852,, 1982. 

Tilling, R., Ridout, A., and Shepher, A.: Estimating Arctic sea ice thickness and volume using CryoSat-2 radar altimeter data, Adv. Space Res., 62, 1203–1225, 2018. 

Toyota, T. and Kimura, N.: An examination of the sea ice rheology for seasonal ice zones based on ice drift and thickness observations, J. Geophys. Res.-Oceans, 123, 1406–1428, 2018. 

Tremblay, L. B. and Mysak, L. A.: Modeling sea ice as a granular material, including the dilatancy effect, J. Phys. Oceanogr., 27, 2342–2360, 1997. 

Tremblay, L. B., and Hakakian, M.: Estimating the sea ice compressive strength from satellite derived sea ice drift and NCEP reanalysis data, J. Phys. Oceanogr., 36, 2165–2172, 2006.  

Tschudi, M., Meier, W., Stewart, J., Fowler, C., and Maslanik, J.: Polar Pathfinder daily 25 km EASE-Grid sea ice motion vectors, version 4, dataset 0116, NASA NationalSnow and Ice Data Center Distributed Active Archive Center, Boulder, CO, USA,, 2019. 

Uotila, P., Farrell, S. O., Marsland, S., and Bi, D.: A sea-ice sensitivity study with a global ocean-ice model, Ocean Model., 51, 1–18,, 2012. 

Vancoppenolle, M., Fichefet, T., Goosse, H., Bouillon, S., Madec, G., and Maqueda, M. A. M.: Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 1. Model description and validation, Ocean Model., 27, 33–53, 2009. 

Wunsch, C.: The Ocean Circulation Inverse Problem, Cambridge Univ. Press, Cambridge, 442 pp., 1996. 

Yaremchuk, M. and Martin, P.: On Sensitivity Analysis within the 4DVAR Framework, Mon. Weather Rev., 142, 774–787, 2014. 

Yaremchuk, M. and Nechaev, D.: Covariance localization with diffusion-based correlation models, Mon. Weather Rev., 141, 848–860, 2013. 

Yaremchuk, M., Nechaev, D., and Panteleev, G.: A method of successive corrections of the control subspace in the reduced-order variational data assimilation, Mon. Weather Rev., 137, 2966–2978, 2009. 

Yaremchuk, M., Townsend, T., Panteleev, G., Hebert, D., and Allard, R.: Advancing short‐term forecasts of ice conditions in the Beaufort Sea, J. Geophys. Res., 124, 807–820,, 2019. 

Zhang J. and Hibler III, W. D.: On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res., 102, 8691–8702, 1997. 

Zhang, J. and Rothrock, D.: Modeling global sea ice with a thickness and enthalpy distribution model in generalized curvilinear coordinates, Mon. Weather Rev., 131, 845–861, 2003. 

Zhang, Y.-F. and Bitz, C. M.: Insights on sea ice data assimilation from perfect model observing system simulation experiments, J. Climate, 31, 5911–5926, 2018. 

Short summary
In the CICE6 community model, rheology and landfast grounding/arching effects are simulated by functions of sea ice thickness and concentration with a set of fixed parameters empirically adjusted to optimize model performance. In this study we consider a spatially variable extension for representing these parameters in the two-dimensional elastic–viscoplastic (EVP) sea ice model and analyze the feasibility of the optimization of these parameters through the 4D-Var data assimilation approach.