Parameter optimization in sea ice models with elastic–viscoplastic rheology

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 k_T 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 (k₁,k₂ 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 $\mathit{\{}{P}^{*},e,k_t,k_{\mathrm{1}},k_{\mathrm{2}}\mathit{\}}$ 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/m² (Hibler and Walsh, 1982), 15–20 kN/m² (Kreyscher et al., 1997, 2000) and 30–45 kN/m² (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 k_T=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 k₁ and maximum basal stress parameter k₂. The optimal values were found to be 8 and 15 N/m³, respectively, with higher sensitivity with respect to k₁ (Lemieux et al., 2015). However, fixed values of k_T,k₁ and k₂ 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:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\rho }}hA({\partial }_t+\stackrel{\mathrm{̃}}{f}\mathit{k}\times )\mathit{u}=\mathrm{div}\mathit{\sigma }+{\mathit{\tau }}_{\mathrm{b}}+{\mathit{\tau }}_{\mathrm{a}}+{\mathit{\tau }}_{\mathrm{w}},{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}{\partial }_t\mathit{\sigma }+e^{\mathrm{2}}\mathit{\sigma }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{tr}\mathit{\sigma }(\mathrm{1}-e^{\mathrm{2}})\mathbf{I}=\\ \text{(2)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{p}}\left[(\mathrm{1}+k_{\mathrm{T}})\dot{\mathit{ϵ}}-{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{2}}}(\mathrm{1}-k_{\mathrm{T}})\mathbf{I}\right]/{\mathrm{\Delta }}^{\ast },\text{(3)}& {\displaystyle }& {\displaystyle }{\partial }_{t}h=\mathrm{div}\left(h\mathit{u}\right),\text{(4)}& {\displaystyle }& {\displaystyle }{\partial }_{t}A=\mathrm{div}\left(A\mathit{u}\right).\end{array}$$

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 $\mathit{u}=\mathit{\{}u,v\mathit{\}}$ are the 2D fields of ice effective thickness, concentration and velocity; and σ and $\dot{\mathit{ϵ}}$ are the 2D fields of ice stress and the deformation rate tensors:

$$\begin{array}{}\text{(5)}& \mathit{\sigma }=\left[\begin{array}{cc}{\mathit{\sigma }}_{xx}& {\mathit{\sigma }}_{xy}\\ {\mathit{\sigma }}_{xy}& {\mathit{\sigma }}_{yy}\end{array}\right];\dot{\mathit{ϵ}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{cc}\mathrm{2}{\partial }_{x}u& {\partial }_{x}v+{\partial }_{y}u\\ {\partial }_{x}v+{\partial }_{y}u& \mathrm{2}{\partial }_{y}v\end{array}\right].\end{array}$$

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

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }\left(\dot{\mathit{ϵ}}\right)={\displaystyle \frac{\mathrm{1}}{e}}{\left[(e^{\mathrm{2}}-\mathrm{1})(\mathrm{tr}\dot{\mathit{ϵ}}{)}^{\mathrm{2}}+\mathrm{2}\mathrm{tr}({\dot{\mathit{ϵ}}}^{\mathrm{2}})\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

To avoid numerical singularities at $\dot{\mathit{ϵ}}=\mathrm{0}$, the values of Δ are limited from below by the additional parameter ${\stackrel{\mathrm{̃}}{\mathrm{\Delta }}}^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ s⁻¹, so that ${\mathrm{\Delta }}^{\ast }\left(\dot{\mathit{ϵ}}\right)=max(\mathrm{\Delta },{\stackrel{\mathrm{̃}}{\mathrm{\Delta }}}^{\ast })$.

The empirical parameters ${\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}},P^{*},k_{\mathrm{T}}$ 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 P_p is related to ice thickness and concentration in accordance with the rheology of Hibler (1979):

$$\begin{array}{}\text{(7)}& P_{\mathrm{p}}=P^{*}hA\exp \left[-\stackrel{\mathrm{̃}}{\mathit{\alpha }}(\mathrm{1}-A)\right].\end{array}$$

The typical values of the ice strength parameter P^* and $\stackrel{\mathrm{̃}}{\mathit{\alpha }}$ 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):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-C_{\mathrm{b}}\mathit{u}\equiv -\mathit{\theta }(Ah-A{h}_{\mathrm{b}}/{\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}){\displaystyle \frac{k_{\mathrm{2}}\mathit{u}}{\left|\mathit{u}\right|+{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0}}}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }\exp \left[-{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{\mathrm{b}}(\mathrm{1}-A)\right],\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{w}}=-{\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{w}}A|\mathit{u}-{\mathit{u}}_{\mathrm{w}}|{\mathbf{R}}_{\mathrm{\Theta }}(\mathit{u}-{\mathit{u}}_{\mathrm{w}}),\end{array}$$

where θ is the Heaviside step function, h_b is the ocean depth, R_Θ is the 2×2 matrix rotating the velocity vector by the turning angle Θ counterclockwise, ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{5}}$ m/s and u_w is the water velocity (set to zero in the present study). The values of other parameters (${\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}},{\mathit{\rho }}_{\mathrm{w}},{\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{\mathrm{b}},{\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ and k₂) are listed in Table 1.

Table 1 Model and assimilation system configuration parameters.

In contrast to the previous studies, where the free empirical parameters $P^{*},e,k_{\mathrm{T}}$ and k₂ were assumed to be constant, the present study attempts to retrieve their spatial variability from synthetic (satellite) observations of the sea ice state vector $\mathit{C}\equiv \mathit{\{}h,A,\mathit{u}\mathit{\}}$ 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 ${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{xx}+{\mathit{\sigma }}_{yy}$, ${\mathit{\sigma }}_{\mathrm{2}}={\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{yy}$, σ₃=σ_xy, ${\dot{\mathit{ϵ}}}_{\mathrm{1}}={\partial }_{x}u+{\partial }_{y}v$, ${\dot{\mathit{ϵ}}}_{\mathrm{2}}={\partial }_{x}u-{\partial }_{y}v$, ${\dot{\mathit{ϵ}}}_{\mathrm{3}}={\partial }_{y}u+{\partial }_{x}v$, bulk viscosity $\mathit{\zeta }=P_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})/\mathrm{2}{\mathrm{\Delta }}^{\ast }$, and replacement pressure $P=P_{\mathrm{p}}\mathrm{\Delta }/{\mathrm{\Delta }}^{\ast }$, Eq. (2) can be split into three decoupled relationships:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\partial }_t{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{1}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}=\mathrm{2}\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{1}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}-P(\mathrm{1}-k_{\mathrm{T}})/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}},\text{(11)}& {\displaystyle }{\partial }_t{\mathit{\sigma }}_{\mathrm{2}}+e^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{2}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}=\mathrm{2}\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{2}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}},\text{(12)}& {\displaystyle }{\partial }_t{\mathit{\sigma }}_{\mathrm{3}}+e^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{3}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}=\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{3}}/{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}.\end{array}$$

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 $\mathit{\delta }{t}_{\mathrm{s}}=\mathit{\epsilon }{\stackrel{\mathrm{̃}}{T}}_{\mathrm{d}}$ (ε=0.00347; see Table 1):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{1}}^s-(\mathrm{1}+\mathit{\epsilon }{)}^{-\mathrm{1}}[{\mathit{\sigma }}_{\mathrm{1}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{1}}/{\mathrm{\Delta }}^{*}\\ \text{(13)}& {\displaystyle }& {\displaystyle }-\mathit{\epsilon }P(\mathrm{1}-k_{\mathrm{T}})]=\mathrm{0},\text{(14)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{2}}^s-(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{1}}[{\mathit{\sigma }}_{\mathrm{2}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{2}}/{\mathrm{\Delta }}^{*}]=\mathrm{0},\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{3}}^s-(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{1}}[{\mathit{\sigma }}_{\mathrm{3}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{3}}/\mathrm{2}{\mathrm{\Delta }}^{*}]=\mathrm{0}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }m{u}^s/\mathit{\delta }{t}_{\mathrm{s}}-\mathrm{mf}{v}^s+C_{\mathrm{b}}{u}^s+{\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{w}}A\left|\mathit{u}\right|(u^s\cos \mathrm{\Theta }-v^s\sin \mathrm{\Theta })\\ \text{(16)}& {\displaystyle }& {\displaystyle }-mu/\mathit{\delta }{t}^s-{\partial }_x{\mathit{\sigma }}_{xx}^s-{\partial }_y{\mathit{\sigma }}_{\mathrm{3}}^s-{\mathit{\tau }}_{\mathrm{a}x}=\mathrm{0},{\displaystyle }& {\displaystyle }m{v}^s/\mathit{\delta }{t}_{\mathrm{s}}+\mathrm{mf}{u}^s+C_{\mathrm{b}}{v}^s+{\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{w}}A\left|\mathit{u}\right|(u^s\sin \mathrm{\Theta }+v^s\cos \mathrm{\Theta })\\ \text{(17)}& {\displaystyle }& {\displaystyle }-mv/\mathit{\delta }{t}^s-{\partial }_y{\mathit{\sigma }}_{yy}^s-{\partial }_x{\mathit{\sigma }}_{\mathrm{3}}^s-{\mathit{\tau }}_{\mathrm{a}y}=\mathrm{0}.\end{array}$$

Here $m=\stackrel{\mathrm{̃}}{\mathit{\rho }}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 u_n and a simplified Lax–Wendroff scheme:

$$\begin{array}{}\text{(18)}& {\displaystyle }{h}^n& {\displaystyle }-h+\mathit{\delta }t\left[\mathrm{div}\right({\mathit{u}}_{n}h)+\mathit{\nu }\widehat{D}h]=\mathrm{0},\text{(19)}& {\displaystyle }{A}^n& {\displaystyle }-A+\mathit{\delta }t\left[\mathrm{div}\right({\mathit{u}}_{n}A)+\mathit{\nu }\widehat{D}A]=\mathrm{0},\end{array}$$

where $\widehat{D}$ is the Laplacian operator and $\mathit{\nu }=\mathit{\delta }t{\mathit{u}}_n^{\mathrm{2}}/\mathrm{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. 13–15) 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

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\mathcal{J}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{\mathrm{\Omega }}\left[W_h(h-h^{\prime }{)}^{\mathrm{2}}+W_A(A-A^{\prime }{)}^{\mathrm{2}}+W_u(\mathit{u}-{\mathit{u}}^{\prime }{)}^{\mathrm{2}}\right.\\ & \left.+{\stackrel{\mathrm{̃}}{W}}_h(\widehat{D}h{)}^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{W}}_A(\widehat{D}A{)}^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{W}}_u(\widehat{D}\mathit{u}{)}^{\mathrm{2}}\right].\end{array}\end{array}$$

Here, W and $\stackrel{\mathrm{̃}}{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 $\mathit{C}=[{\mathit{C}}_{\mathrm{0}},{\mathit{C}}_{\mathrm{p}}]$, which includes the initial state of the model ${\mathit{C}}_{\mathrm{0}}\equiv \mathit{X}{|}_{t=\mathrm{0}}$, and other control fields C_p, 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 C_p 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 $\widehat{\mathit{X}}$ (adjoints of the state variables X) in Ω (e.g., Le Dimet and Talagrand, 1986) and minimize the modified cost function

$$\begin{array}{}\text{(21)}& {\mathcal{J}}_m=\mathcal{J}+\widehat{\mathcal{J}}\equiv \mathcal{J}\left(\mathit{X}\right)+\mathit{M}(\mathit{X},\mathit{C}{)}^{\mathsf{T}}\widehat{\mathit{X}},\end{array}$$

with respect to $\mathit{X},\widehat{\mathit{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 $\mathit{\delta }\mathit{X},\mathit{\delta }\widehat{\mathit{X}}$ and δC) to zero:

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }\widehat{\mathit{X}}}}=\mathit{M}(\mathit{X},\mathit{C})=\mathrm{0},\text{(23)}& {\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }\mathit{X}}}={\displaystyle \frac{\mathit{\delta }\mathcal{J}}{\mathit{\delta }\mathit{X}}}+{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }\mathit{X}}}({\mathit{M}}_X\mathit{\delta }\mathit{X}{)}^{\mathsf{T}}\widehat{\mathit{X}}={\displaystyle \frac{\mathit{\delta }\mathcal{J}}{\mathit{\delta }\mathit{X}}}+{\mathit{M}}_X^{\mathsf{T}}\widehat{\mathit{X}}=\mathrm{0},\text{(24)}& {\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }\mathit{C}}}={\displaystyle \frac{\mathit{\delta }\widehat{\mathcal{J}}}{\mathit{\delta }\mathit{M}}}{\displaystyle \frac{\mathit{\delta }\mathit{M}}{\mathit{\delta }\mathit{C}}}=\mathrm{0}.\end{array}$$

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 $\widehat{\mathcal{J}}$ with respect to X and involves the product of the model operator M_X 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 M_X) forced by the derivatives of 𝒥 to obtain the values of the adjoint variables $\widehat{\mathit{X}}$. Finally, the derivatives of 𝒥_m are computed from Eq. (24) using the chain rule and the values of X and $\widehat{\mathit{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 $\mathit{C}=\mathit{\left\{}\mathit{u}\right(\mathrm{0},\mathit{x}),A(\mathrm{0},\mathit{x}),h(\mathrm{0},\mathit{x}),{\mathit{C}}_{\mathrm{p}}(t,\mathit{x}\left)\mathit{\right\}}$ at the nth iteration;

2. run the forward model (Eqs. 13–19) to compute X_n, the derivatives δ𝒥∕δX and the value of 𝒥_n (Eq. 20);

3. run the adjoint model forced by $-\mathit{\delta }\mathcal{J}/\mathit{\delta }\mathit{X}$ (cf. Eq. 23) to obtain ${\widehat{\mathit{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. 13–19), the outlined optimization requires development of the routines for computing the gradients as well as the tangent linear model M_X and its adjoint ${\mathit{M}}_X^{\mathsf{T}}$. 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: http://autodiff.com/tamc, last access: 1 December 2020; OpenAD: https://www.mcs.anl.gov/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 M_X 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 M_X 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 M_X 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 ${\widehat{\mathit{X}}}^{\mathsf{T}}{\mathit{M}}_X\mathit{\delta }\mathit{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 δt_s 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 δt_s 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 δt_s 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 δt_s 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 M_X.

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 δt_s 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 ${\mathit{\epsilon }}_{\mathrm{N}}{\mathit{\sigma }}_i,i=\mathrm{1},\mathrm{2},\mathrm{3}$ appear inside the square brackets of the linearized Eqs. (13)–(15), where ε_N is the ratio of δt_s to the Newtonian damping timescale T_N (see Appendix A for details). Numerical experimentation has shown that this approach worked generally well using the Newtonian damping timescale T_N of 7δt_s.

Additional experiments have shown that T_N must decrease to 3δt_s 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 $\mathit{C}=\mathit{\left\{}h\right(\mathit{x},\mathrm{0}),A(\mathit{x},\mathrm{0}\left)\mathit{\right\}}$ were slightly perturbed by the realizations of a random function R (viz. $\mathit{C}(\mathit{x},\mathrm{0})\to \mathit{C}(\mathit{x},\mathrm{0})+\mathit{ϵ}R\left(\mathit{x}\right)\equiv \mathit{C}+\mathit{\delta }\mathit{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:

$$\begin{array}{}\text{(25)}& \mathrm{\Phi }\left(\mathit{ϵ}\right)={\displaystyle \frac{\left|\mathit{X}\right(\mathit{C}+\mathit{\delta }\mathit{C})-\mathit{X}(\mathit{C})-TL(\mathit{\delta }\mathit{C}\left)\right|}{\left|\mathit{X}\right(\mathit{C}\left)\right|}}\end{array}$$

, using the initial conditions listed in the argument, and $|\cdot |$ 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 1 TL 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 ϵ² 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, k_T and k₂. 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 (http://www.cpom.ucl.ac.uk/csopr/seaice.html, 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 (https://icesat-2.gsfc.nasa.gov/, 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 ($\mathrm{0}\le A\le \mathrm{1}$) and the control fields of the rheological parameters (listed at the bottom of the left column in Table 1).

Table 2 List of the performed experiments

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 k₂ and k_T. The control vector in these experiments included only parameters k₂ and k_T, 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 (δx_p=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 $\mathrm{36}=(\mathrm{75}/\mathrm{15}+\mathrm{1})(\mathrm{30}/\mathrm{15}+\mathrm{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 k_T

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 k_T 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 2 Results of the OSSEs optimizing k_T. (a, b) True SIC and SIV with k_T=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 k_T only (solid line) and with the joint optimization of k_T, h and A (dashed line) (e). Left part of panel (f) shows optimized k_T for the experiment with true first-guess SIC/SIT distributions at t=0. (g, h) Results of the OSSEs optimizing k_T with true k_T=0.2 and zonal wind of 6 m/s: optimized velocity field and concentration at t=3 d (g); optimized k_T (h).

The first-guess solution was forced by the same wind but with k_T=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 (k_T=0).

The 4D-Var optimization of only k_T (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 k_T (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 k_T enabled formation of landfast ice in the western part of the model domain.

Figure 3 Results of the k₂ optimization with ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}=\mathrm{0.8}$. (a, b) True SIC and SIV with k₂ = 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 k₂ = 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 k₂ (f) and the optimized k₂ distribution (g).

Obviously, a similar effect could be achieved with much higher values of k_T. To remove this ambiguity, we added an additional term into the cost function which is proportional to the integral of $k_{\mathrm{T}}^{\mathrm{2}}$ over the domain. By minimizing the magnitude of k_T, we find the minimum value of k_T necessary for holding ice in place. Note that optimization of k_T 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 k_T. 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 k_T specified over the entire Arctic Ocean and found that k_T 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 k_T 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 k_T=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 k_T 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 k_T values is very close to the true value of k_T = 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 k_T only in the key regions, and thus there is no need to specify uniform k_T 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 k₂

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 k_T were set to zero, while the following true values of the RPs (Lemieux et al., 2016) were used: the critical thickness parameter ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ = 8, the basal stress parameter ${\stackrel{\mathrm{̃}}{\mathit{\alpha }}}_{\mathrm{b}}$ = 20 and the maximum basal stress parameter k₂ = 15 N/m³. 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 k₂ 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 k₂, 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 k₂ is shown in Fig. 3g. The maximum values of k₂ are very close to the true k₂ = 15 N/m³. Note that the true k₂ value was specified ad hoc and the grounding effect formally could be reached with smaller values of k₂. This is clearly demonstrated in Fig. 3g, where k₂ is about 12 N/m³ 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 k₂ was added to the cost function, and the optimized field of k₂ 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 k₂, which should invoke faster convergence of the minimization procedure.

To demonstrate the robustness of the k₂ optimization with respect to possible variations in ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$, we conducted a similar experiment with a smaller true ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ = 2.5. As seen in Fig. 4a and b, the decrease in ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{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 k₂ = 0, i.e., all ice moving eastward. However, after 4D-Var DA and optimization of k₂, the reconstructed solution (Fig. 4d) is nearly identical to the true solution (Figure 4b). This optimized k₂ (Fig. 4c) has the same spatial structure but a slightly smaller (∼  14 N/m²) maximum value compared to the experiment where ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ = 8.

Figure 4 Results of k₂ optimization with ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}=\mathrm{2.5}$. (a, b) True SIC (color) and SIV distributions with k₂=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 k₂ distribution (c) and optimized SIC and SIV at t = 3 d (d).

The parameterization of grounding landfast ice also includes the critical thickness parameter ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{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 ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ than k₂ (Lemieux et al., 2015) because ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}$ can be interpreted as a scaling coefficient in the definition of the critical ice thickness ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}{h}_{\mathrm{c}}=A{h}_{\mathrm{b}}$ (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 k₁ 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 k_T and k₂.

Figure 5 True 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 P_tr 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 ², 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/m${}^{\mathrm{2}}\le P^{*}\le \mathrm{32.5}$ kN/m² and $\mathrm{1}\le e\le$ 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/m². 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 ${\mathit{\sigma }}_I\equiv P_{\mathrm{tr}}=-\mathrm{tr}\mathit{\sigma }/\mathrm{2}$ at the beginning of the true state integration and after 3 d. The P_tr distribution has a clear maximum near the location with the coordinates (500 km, 500 km), which corresponds to the maximum pressure (∼ 40 kN/m²) 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, P_tr 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/m² 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/m², e = 2 and true wind forcing.

Figure 6 Results of the GYRE-0 OSSE. Suboptimal distribution of SIT, SIC (a) and P_tr0 (c) for t = 3 d after optimization of the initial conditions u₀, h₀ and A₀ using the first-guess values of P^* = 27 kN/m² and e = 2. The deviation norms S_u and S_tr 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/m²) 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 ${\mathit{C}}_{\mathrm{0}}=\mathit{\{}{\mathit{u}}_{\mathrm{0}},h_{\mathrm{0}},A_{\mathrm{0}}\mathit{\}}$. Then we sequentially optimized rheological components of the control vector ${\mathit{C}}_{\mathrm{rh}}=\mathit{\{}{P}^{*},e\mathit{\}}$ and finally conducted an additional optimization of the full control vector $\mathit{C}=\mathit{\{}{\mathit{C}}_{\mathrm{0}},{\mathit{C}}_{\mathrm{rh}}\mathit{\}}$. 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 C₀ with fixed P^* and e (left panels) and after additional spatial optimization of the P^* and e (right panels). The major result of optimizing C_rh is an improvement of the SIV fields. The formal quantitative measure of the SIV improvement was evaluated by the function

$${\displaystyle }{S}_{\mathrm{u}}=\sum _{\mathrm{\Omega },k}|{\mathit{u}}_{\mathrm{opt}}^k-{\mathit{u}}_{\mathrm{true}}|,$$

where Ω is the model domain and index k = 1,2 enumerates suboptimal optimization stages: k=1 using the initial conditions control vector C₀ only and k=2 employing the full control vector $\mathit{C}=\mathit{\{}{\mathit{C}}_{\mathrm{0}},{\mathit{C}}_{\mathrm{rh}}\mathit{\}}$. It was found that S_u reduced almost 1.5 times after the additional optimization of the rheological parameters C_rh (Fig. 6a, b). Visual comparison of the suboptimal and fully optimized SIC shows a certain improvement after C_rh 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 7 Results of the GYRE-W OSSE: optimized distribution of SIT, SIC (a) and P_tr (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/m², left) and e.

In this experiment, we found that optimization of C_rh yields only a marginal improvement of the SIT distribution. In particular, SD$(h_{\mathrm{opt}}^k-h_{\mathrm{true}})$ decreased from 0.23 m, after optimization of the initial conditions, to 0.2 m, after additional C_rh optimization. The minor impact of C_rh 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 C_rh over the relatively short timescale of the 3 d assimilation window.

As expected, optimization of the rheological parameters C_rh 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 P_tr decreased by ∼ 35 %, after the additional RP optimization. Note, also, that the fully optimized P_tr maximum (Fig. 6d) is located in close proximity to the true P_tr maximum (Fig. 5d), while the maximum in the suboptimal P_tr 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 P_tr 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 $\mathit{C}=\mathit{\{}{\mathit{C}}_{\mathrm{0}},{\mathit{C}}_{\mathrm{rh}}\mathit{\}}$ 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/m², 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 S_u increased five times up to 0.64 m/s.

At the same time, degradation of the SIT retrieval was not as significant, with SD$(h_{\mathrm{opt}}^k-h_{\mathrm{true}})$ 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 P_tr 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 P_tr 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.

4.2 PIZ OSSE

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 A₀ from the control vector C₀. 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/m². 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 8 True 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 P_tr of the internal stress tensor. White arrows show the wind stress forcing at t=0. The distributions of P^* (kN/m²) 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 P_tr, shown in the lower panels of Fig. 8, has two clear maxima in the south, with the magnitudes of about 70 and 50 kN/m². Note that due to ice convergence causing SIT growth both P_tr 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/m², e = 2 and the exact wind forcing. The optimization was conducted in three steps: by first optimizing C₀, then optimizing C_rh, and, finally, simultaneously optimizing C₀ and C_rh.

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 9 Results of the PIZ OSSE: suboptimal distributions of SIT (a) and P_tr (c) at t = 3 d after optimizing the initial conditions ${\mathit{u}}_{\mathrm{0}},h_{\mathrm{0}},A_{\mathrm{0}}$ and the first-guess values of P^* = 27.5 kN/m² and e = 2. The discrepancies of S_u and S_tr 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/m²) and e. Panels (g) and (f) show optimized distribution of P^* (kN/m²) and e for the experiment with moderate wind of 10 m/s.

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

Additional optimization of the rheological parameters, C_rh, significantly improves the reconstructed SIV practically everywhere, with the formal measure of the uncertainty S_u decreasing by almost one-half from 56 to 30 m/s (Fig. 9b). Similar improvements are visible in the P_tr distribution (Fig. 9d). The SD[P_tr(opt)−P_tr(true)] decreased to 7.4 kN/m² (by 40 %), and the fully optimized P_tr has two maxima as in the true P_tr 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/m². 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 k₁), 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 k₂ and k_T, 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 k_T could be achieved by specifying k_T 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 k_T can be reduced to specifying nonzero k_T only within these localized regions of the domain. Interestingly, the optimization of both k₂ and k₁ 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 k₁, which confines k₂ variability to shallow areas and narrow straits areas in a high-resolution pan-Arctic sea ice model. It should be noted that retrieval of k₁ is not straightforward because of nondifferentiability of the grounding parameterization scheme (Eq. 8). Optimizing k₁ 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 k₁, 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 P_tr fields. Accurate forecasting of P_tr 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 P_tr. OSSEs with weaker winds (∼ 0.0–0.15 N/m²) 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 k₂ and k_T remained virtually unchanged. The incoming satellite platforms (e.g., ICESat-2, https://icesat-2.gsfc.nasa.gov, 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, http://autodiff.com/tamc/, last access: 1 December 2020; OpenAD, https://www.mcs.anl.gov/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 ${\dot{\mathit{ϵ}}}_{\mathrm{1},\mathrm{2},\mathrm{3}},m,\mathit{\zeta },P_{\mathrm{p}}$, P and C_b of X are given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{1}}={\partial }_x\mathit{\delta }u+{\partial }_y\mathit{\delta }v;\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{2}}={\partial }_x\mathit{\delta }u-{\partial }_y\mathit{\delta }v;\\ \text{(A1)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{3}}={\partial }_x\mathit{\delta }v+{\partial }_y\mathit{\delta }u;\mathit{\delta }m=\stackrel{\mathrm{̃}}{\mathit{\rho }}(A\mathit{\delta }h+h\mathit{\delta }A),{\displaystyle }& {\displaystyle }\mathit{\delta }\mathit{\zeta }=\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}}(\mathrm{1}+k_{\mathrm{T}})/\mathrm{2}{\mathrm{\Delta }}^{\ast }+\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{T}}}{P}_{\mathrm{p}}/\mathrm{2}{\mathrm{\Delta }}^{\ast }\\ \text{(A2)}& {\displaystyle }& {\displaystyle }-\underset{\mathrm{¯}}{\mathit{\delta }{\mathrm{\Delta }}^{\ast }}{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})/\mathrm{2}{\mathrm{\Delta }}^{*\mathrm{2}},{\displaystyle }& {\displaystyle }\mathit{\delta }{P}_{\mathrm{p}}=(P_{\mathrm{p}}^{*}h\mathit{\delta }A+P_{\mathrm{p}}^{*}A\mathit{\delta }h+Ah\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}^{*}})\exp (-\stackrel{\mathrm{̃}}{\mathit{\alpha }}(\mathrm{1}-A\left)\right)\\ \text{(A3)}& {\displaystyle }& {\displaystyle }+P_{\mathrm{p}}\stackrel{\mathrm{̃}}{\mathit{\alpha }}\mathit{\delta }A,\text{(A4)}& {\displaystyle }& {\displaystyle }\mathit{\delta }P={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Delta }}^{\ast }}}\left[\mathrm{\Delta }\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}}+P_{\mathrm{p}}\mathit{\delta }\mathrm{\Delta }-P_{\mathrm{p}}{\displaystyle \frac{\mathrm{\Delta }}{{\mathrm{\Delta }}^{\ast }}}\underset{\mathrm{¯}}{\mathit{\delta }{\mathrm{\Delta }}^{\ast }}\right],\text{(A5)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{\mathrm{\Delta }}^{\ast }=\mathit{\theta }(\mathrm{\Delta }-{\stackrel{\mathrm{̃}}{\mathrm{\Delta }}}^{\ast })\underset{\mathrm{¯}}{\mathit{\delta }\mathrm{\Delta }},\text{(A6)}& {\displaystyle }& {\displaystyle }\mathit{\delta }\mathrm{\Delta }={\displaystyle \frac{\mathrm{1}}{e^{\mathrm{2}}\mathrm{\Delta }}}\left[e^{\mathrm{2}}{\dot{\mathit{ϵ}}}_{\mathrm{1}}\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{1}}+{\dot{\mathit{ϵ}}}_{\mathrm{2}}\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{3}}-{\displaystyle \frac{\mathit{\delta }e}{e}}({\dot{\mathit{ϵ}}}_{\mathrm{2}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}^{\mathrm{2}})\right],\text{(A7)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{C}_{\mathrm{b}}=C_{\mathrm{b}}\left\{{\displaystyle \frac{\mathit{\delta }h}{h_{\mathrm{b}}}}-{\displaystyle \frac{\mathit{\delta }A}{{\stackrel{\mathrm{̃}}{k}}_{\mathrm{1}}}}+{\displaystyle \frac{\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{2}}}}{k_{\mathrm{2}}}}-{\displaystyle \frac{\mathit{u}\cdot \mathit{\delta }\mathit{u}}{{\mathit{u}}^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{u}}_{\mathrm{0}}\left|\mathit{u}\right|}}\right\}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\delta }{\mathit{\sigma }}_{\mathrm{1}}^s=(\mathrm{1}+\mathit{\epsilon }{)}^{-\mathrm{1}}[\mathit{\delta }{\mathit{\sigma }}_{\mathrm{1}}-\mathit{\epsilon }\underset{\mathrm{¯}}{\mathit{\delta }P}(\mathrm{1}-k_{\mathrm{T}})\\ {\displaystyle }& {\displaystyle }+\mathit{\epsilon }P\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{T}}}+\mathit{\epsilon }\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{1}}/{\mathrm{\Delta }}^{\ast }+\mathit{\epsilon }{P}_{\mathrm{p}}{\dot{\mathit{ϵ}}}_{\mathrm{1}}\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{T}}}/{\mathrm{\Delta }}^{\ast }+\\ {\displaystyle }& {\displaystyle }+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{1}}/{\mathrm{\Delta }}^{\ast }-\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{1}}\underset{\mathrm{¯}}{\mathit{\delta }{\mathrm{\Delta }}^{\ast }}/{\mathrm{\Delta }}^{\ast \mathrm{2}}\\ \text{(A8)}& {\displaystyle }& {\displaystyle }-\mathit{\epsilon }{}_{\mathrm{N}}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{1}}],{\displaystyle }& {\displaystyle }\mathit{\delta }{\mathit{\sigma }}_{\mathrm{2}}^s=(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{1}}[\mathit{\delta }{\mathit{\sigma }}_{\mathrm{2}}+\mathit{\epsilon }\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{2}}/{\mathrm{\Delta }}^{\ast }\\ {\displaystyle }& {\displaystyle }+\mathit{\epsilon }{P}_{\mathrm{p}}{\dot{\mathit{ϵ}}}_{\mathrm{2}}\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{T}}}/{\mathrm{\Delta }}^{\ast }+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{2}}/{\mathrm{\Delta }}^{\ast }-\\ {\displaystyle }& {\displaystyle }-\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{2}}\underset{\mathrm{¯}}{\mathit{\delta }{\mathrm{\Delta }}^{\ast }}/{\mathrm{\Delta }}^{\ast \mathrm{2}}-\mathit{\epsilon }{}_{\mathrm{N}}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{2}}],\\ \text{(A9)}& {\displaystyle }& {\displaystyle }-\mathrm{2}\mathit{\epsilon }e\underset{\mathrm{¯}}{\mathit{\delta }e}(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{2}}({\mathit{\sigma }}_{\mathrm{2}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{2}}/{\mathrm{\Delta }}^{\ast }){\displaystyle }& {\displaystyle }\mathit{\delta }{\mathit{\sigma }}_{\mathrm{3}}^s=(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{1}}[\mathit{\delta }{\mathit{\sigma }}_{\mathrm{3}}+\mathit{\epsilon }\underset{\mathrm{¯}}{\mathit{\delta }{P}_{\mathrm{p}}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{3}}/\mathrm{2}{\mathrm{\Delta }}^{\ast }\\ {\displaystyle }& {\displaystyle }+\mathit{\epsilon }{P}_{\mathrm{p}}\underset{\mathrm{¯}}{\mathit{\delta }{k}_{\mathrm{T}}}{\dot{\mathit{ϵ}}}_{\mathrm{3}}/\mathrm{2}{\mathrm{\Delta }}^{\ast }+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})\mathit{\delta }{\dot{\mathit{ϵ}}}_{\mathrm{3}}/\mathrm{2}{\mathrm{\Delta }}^{\ast }-\\ {\displaystyle }& {\displaystyle }-\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{3}}\underset{\mathrm{¯}}{\mathit{\delta }{\mathrm{\Delta }}^{\ast }}/\mathrm{2}{\mathrm{\Delta }}^{\ast \mathrm{2}}-\mathit{\epsilon }{}_{\mathrm{N}}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{3}}]\\ \text{(A10)}& {\displaystyle }& {\displaystyle }-\mathrm{2}\mathit{\epsilon }e\underset{\mathrm{¯}}{\mathit{\delta }e}(\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }{)}^{-\mathrm{2}}({\mathit{\sigma }}_{\mathrm{3}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}}){\dot{\mathit{ϵ}}}_{\mathrm{3}}/\mathrm{2}{\mathrm{\Delta }}^{\ast }),{\displaystyle }& {\displaystyle }(m\mathit{\delta }{u}^s+u^s\mathit{\delta }m)/\mathit{\delta }{t}_{\mathrm{s}}-f(v^s\mathit{\delta }m+m\mathit{\delta }{v}^s)+\underset{\mathrm{¯}}{\mathit{\delta }{C}_{\mathrm{b}}}{u}^s\\ {\displaystyle }& {\displaystyle }+C_{\mathrm{b}}\mathit{\delta }{u}^s+{\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{w}}A\left|\mathit{u}\right|(\mathit{\delta }{u}^s\cos \mathrm{\Theta }-\mathit{\delta }{v}^s\sin \mathrm{\Theta })=\\ {\displaystyle }& {\displaystyle }={\mathit{\tau }}_{wx}(u\mathit{\delta }u/{\mathit{u}}^{\mathrm{2}}+\mathit{\delta }A/A)+(u\mathit{\delta }m+m\mathit{\delta }u)/\mathit{\delta }{t}_{\mathrm{s}}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{¯}}{\mathit{\delta }{\mathit{\tau }}_{\mathrm{a}x}}+{\partial }_x\mathit{\delta }{\mathit{\sigma }}_{xx}^s+{\partial }_y\mathit{\delta }{\mathit{\sigma }}_{\mathrm{3}}^s,{\displaystyle }& {\displaystyle }(m\mathit{\delta }{v}^s+v^s\mathit{\delta }m)/\mathit{\delta }{t}_{\mathrm{s}}+f(u^s\mathit{\delta }m+m\mathit{\delta }{u}^s)+\underset{\mathrm{¯}}{\mathit{\delta }{C}_{\mathrm{b}}}{v}^s\\ {\displaystyle }& {\displaystyle }+C_{\mathrm{b}}\mathit{\delta }{v}^s+{\stackrel{\mathrm{̃}}{C}}_{\mathrm{w}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{w}}A\left|\mathit{u}\right|(\mathit{\delta }{u}^s\sin \mathrm{\Theta }+\mathit{\delta }{v}^s\cos \mathrm{\Theta })=\\ {\displaystyle }& {\displaystyle }={\mathit{\tau }}_{wy}(v\mathit{\delta }v/{\mathit{u}}^{\mathrm{2}}+\mathit{\delta }A/A)+(v\mathit{\delta }m+m\mathit{\delta }v)/\mathit{\delta }{t}_{\mathrm{s}}\\ \text{(A12)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{¯}}{\mathit{\delta }{\mathit{\tau }}_{\mathrm{a}y}}+{\partial }_y\mathit{\delta }{\mathit{\sigma }}_{yy}^s+{\partial }_x\mathit{\delta }{\mathit{\sigma }}_{\mathrm{3}}^s,{\displaystyle }& {\displaystyle }\mathit{\delta }{h}^n=\mathit{\delta }t\mathrm{div}(\mathit{\delta }{\mathit{u}}^{n}h+{\mathit{u}}^n\mathit{\delta }h)\\ \text{(A13)}& {\displaystyle }& {\displaystyle }+\mathit{\delta }{t}^{\mathrm{2}}({\mathit{u}}^{\mathrm{2}}\widehat{D}\mathit{\delta }h+\mathit{\delta }({\mathit{u}}^{\mathrm{2}}\left)\widehat{D}h\right)/\mathrm{2},{\displaystyle }& {\displaystyle }\mathit{\delta }{A}^n=\mathit{\delta }t\mathrm{div}(\mathit{\delta }{\mathit{u}}^{n}A+{\mathit{u}}^n\mathit{\delta }A)\\ \text{(A14)}& {\displaystyle }& {\displaystyle }+\mathit{\delta }{t}^{\mathrm{2}}({\mathit{u}}^{\mathrm{2}}\widehat{D}\mathit{\delta }A+\mathit{\delta }({\mathit{u}}^{\mathrm{2}}\left)\widehat{D}A\right)/\mathrm{2}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }{P}_{\mathrm{p}}^{*}}}=-\mathit{\epsilon }{\mathcal{I}}^{\mathsf{T}}\sum _t(\mathrm{1}+k_{\mathrm{T}})Ah\exp (-\stackrel{\mathrm{̃}}{\mathit{\alpha }}(\mathrm{1}-A\left)\right)\\ \text{(A15)}& {\displaystyle }& {\displaystyle }({\dot{\mathit{ϵ}}}_{\mathrm{1}}{\widehat{\mathit{\sigma }}}_{\mathrm{1}}+{\dot{\mathit{ϵ}}}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}{\widehat{\mathit{\sigma }}}_{\mathrm{3}}/\mathrm{2})/{\mathrm{\Delta }}^{\ast },{\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }e}}=\mathit{\epsilon }{\mathcal{I}}^{\mathsf{T}}\sum _t\left\{{\displaystyle \frac{\mathrm{2}e}{(\mathrm{1}+\mathit{\epsilon }{e}^{\mathrm{2}}{)}^{\mathrm{2}}}}[{\mathit{\sigma }}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{3}}{\widehat{\mathit{\sigma }}}_{\mathrm{3}}+\mathit{\epsilon }{P}_{\mathrm{p}}(\mathrm{1}+k_{\mathrm{T}})\right.\\ {\displaystyle }& {\displaystyle }({\dot{\mathit{ϵ}}}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}{\widehat{\mathit{\sigma }}}_{\mathrm{3}}/\mathrm{2})/{\mathrm{\Delta }}^{\ast }]-{\displaystyle \frac{P_{\mathrm{p}}(\mathrm{1}-k_{\mathrm{T}})}{e^{\mathrm{3}}\mathrm{\Delta }(\mathrm{1}+\mathit{\epsilon })}}({\dot{\mathit{ϵ}}}_{\mathrm{2}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}^{\mathrm{2}}){\widehat{\mathit{\sigma }}}_{\mathrm{1}}+\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{P_{\mathrm{p}}}{e^{\mathrm{3}}{\mathrm{\Delta }}^{\ast \mathrm{2}}}}\mathit{\theta }(\mathrm{\Delta }-{\stackrel{\mathrm{̃}}{\mathrm{\Delta }}}^{\ast })({\dot{\mathit{ϵ}}}_{\mathrm{2}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{3}}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\left[(\mathrm{1}+k_{\mathrm{T}})\left({\displaystyle \frac{{\dot{\mathit{ϵ}}}_{\mathrm{1}}{\widehat{\mathit{\sigma }}}_{\mathrm{1}}}{\mathrm{1}+\mathit{\epsilon }}}+{\displaystyle \frac{{\dot{\mathit{ϵ}}}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{3}}/\mathrm{2}}{\mathrm{1}+e^{\mathrm{2}}\mathit{\epsilon }}}\right)\right.\\ \text{(A16)}& {\displaystyle }& {\displaystyle }\left.\left.+(\mathrm{1}-k_{\mathrm{T}}){\displaystyle \frac{{\mathrm{\Delta }}^{\ast }{\widehat{\mathit{\sigma }}}_{\mathrm{1}}}{\mathrm{1}+\mathit{\epsilon }}}\right]\right\},{\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }{k}_{\mathrm{T}}}}=-\mathit{\epsilon }{\mathcal{I}}^{\mathsf{T}}\sum _t{P}_{\mathrm{p}}\left[\right({\dot{\mathit{ϵ}}}_{\mathrm{1}}+\mathrm{\Delta }){\widehat{\mathit{\sigma }}}_{\mathrm{1}}+{\dot{\mathit{ϵ}}}_{\mathrm{2}}{\widehat{\mathit{\sigma }}}_{\mathrm{2}}\\ \text{(A17)}& {\displaystyle }& {\displaystyle }+{\dot{\mathit{ϵ}}}_{\mathrm{3}}{\widehat{\mathit{\sigma }}}_{\mathrm{3}}/\mathrm{2}]/{\mathrm{\Delta }}^{\ast },\text{(A18)}& {\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }{k}_{\mathrm{2}}}}={\mathcal{I}}^{\mathsf{T}}\sum _t{C}_{\mathrm{b}}(\mathit{u}\cdot \widehat{\mathit{u}})/k_{\mathrm{2}},\text{(A19)}& {\displaystyle }& {\displaystyle \frac{\mathit{\delta }{\mathcal{J}}_m}{\mathit{\delta }{\mathit{\tau }}_{\mathrm{a}}}}=-{\mathcal{I}}^{\mathsf{T}}\widehat{\mathit{u}}.\end{array}$$

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 $\mathit{\epsilon }{}_{\mathrm{N}}\mathit{\delta }{\mathit{\sigma }}_{\mathrm{1},\mathrm{2},\mathrm{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. 22–24). 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 M_X 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 σ,

$$\begin{array}{}\text{(B1)}& \mathit{\sigma }={\displaystyle \frac{\mathrm{1}}{e^{\mathrm{2}}}}\left[\mathit{P}-\mathbf{I}{\displaystyle \frac{\mathrm{1}-e^{\mathrm{2}}}{\mathrm{2}}}\mathrm{tr}\mathit{P}\right],\end{array}$$

and substituting (B1) into the momentum equation:

$$\begin{array}{}\text{(B2)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\rho }}hA({\partial }_t+f\mathit{k}\times )\mathit{u}=\mathrm{div}\mathit{\sigma }+{\mathit{\tau }}_{\mathrm{b}}+{\mathit{\tau }}_{\mathrm{a}}+{\mathit{\tau }}_{\mathrm{w}},\text{(B3)}& {\displaystyle }& {\displaystyle }{\partial }_{t}h=\mathrm{div}\left(h\mathit{u}\right),\text{(B4)}& {\displaystyle }& {\displaystyle }{\partial }_{t}A=\mathrm{div}\left(A\mathit{u}\right),\end{array}$$

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

$$\begin{array}{}\text{(B5)}& \mathit{P}=P_{\mathrm{p}}\left[(\mathrm{1}+k_{\mathrm{T}})\dot{\mathit{ϵ}}-{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{2}}}(\mathrm{1}-k_{\mathrm{T}})\mathbf{I}\right]/{\mathrm{\Delta }}^{\ast }.\end{array}$$

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 $\dot{\mathit{ϵ}}$ 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.