Research article 04 Dec 2020
Research article  04 Dec 2020
Parameter optimization in sea ice models with elastic–viscoplastic rheology
 ^{1}Naval Research Laboratory, Stennis Space Center, MS, USA
 ^{2}University of Colorado, Aurora, CO, USA
 ^{3}University of Hawai`i at Mānoa, Honolulu, HI, USA
 ^{1}Naval Research Laboratory, Stennis Space Center, MS, USA
 ^{2}University of Colorado, Aurora, CO, USA
 ^{3}University of Hawai`i at Mānoa, Honolulu, HI, USA
Correspondence: Gleb Panteleev (gleb.panteleev@nrlssc.navy.mil)
Hide author detailsCorrespondence: Gleb Panteleev (gleb.panteleev@nrlssc.navy.mil)
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 twodimensional 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.
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 quasiNewtonian 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 fourdimensional variational (4DVar) 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 3DVar, ensemblebased DA systems utilizing the CICE model (e.g., Zhang and Bitz, 2018). However, to our knowledge, there is no 4DVar 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 socalled 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_{1},k_{2} 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^{2} (Hibler and Walsh, 1982), 15–20 kN/m^{2} (Kreyscher et al., 1997, 2000) and 30–45 kN/m^{2} (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 coarseresolution panArctic CICENEMO 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_{1} and maximum basal stress parameter k_{2}. The optimal values were found to be 8 and 15 N/m^{3}, respectively, with higher sensitivity with respect to k_{1} (Lemieux et al., 2015). However, fixed values of k_{T},k_{1} and k_{2} 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 trialanderror 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 abovelisted algorithms are not suitable for simultaneous optimization of other model parameters such as initial and openboundary 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 shortterm (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 onedimensional (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 shortrange 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 4DVar 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 4DVar 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 twodimensional (2D) sea ice model based on the EVP rheology formulation of Lemieux et al. (2016). Our analysis is based on application of the 4DVar and Observing System Simulation Experiment (OSSE) (e.g., Nitta, 1969; Arnold and Dey, 1986; Nichols 2003, 2010) and follows the conventional twindata 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 syntheticaperture 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 shortterm 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 firstguess 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.
This section provides details of the sea ice model formulation and its associated linearizations, outlines the variational data assimilation system used for optimizing model parameters, and describes synthetic observations used to do so. To distinguish between the parameter fields spatially varying in 2D and the fixed parameter values that were not subject to optimization, the latter are marked by tildes. The major parameters of the model are listed in Table 1.
2.1 EVP sea ice model
2.1.1 Formulation
In the present study, we employed the sea ice model formulation of Lemieux et al. (2016) with the basal stress parameterization and generalized Hibler (1979) yield curve. Equations of the model describe EVP ice physics coupled with sea ice dynamics, which is forced by the stresses τ exerted on ice through its interaction with the bottom τ_{b}, atmosphere τ_{a} and the ocean τ_{w}:
Here div and tr are the divergence and trace operators; k is the vertical unit vector; I is the 2 × 2 identity matrix; h,A and $\mathit{u}=\mathit{\{}u,v\mathit{\}}$ are the 2D fields of ice effective thickness, concentration and velocity; and σ and $\dot{\mathit{\u03f5}}$ are the 2D fields of ice stress and the deformation rate tensors:
The scalar field Δ used for normalizing the righthand side (rhs) in (2) is computed using the following expression (e.g., Hunke, 2001):
To avoid numerical singularities at $\dot{\mathit{\u03f5}}=\mathrm{0}$, the values of Δ are limited from below by the additional parameter ${\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}}^{\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\ast}={\mathrm{10}}^{\mathrm{10}}$ s^{−1}, so that ${\mathrm{\Delta}}^{\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\ast}\left(\dot{\mathit{\u03f5}}\right)=max(\mathrm{\Delta},{\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}}^{\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\ast})$.
The empirical parameters ${\stackrel{\mathrm{\u0303}}{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):
The typical values of the ice strength parameter P^{*} and $\stackrel{\mathrm{\u0303}}{\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):
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{\u0303}}{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{\u0303}}{C}}_{\mathrm{w}},{\mathit{\rho}}_{\mathrm{w}},{\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}}_{\mathrm{b}},{\stackrel{\mathrm{\u0303}}{k}}_{\mathrm{1}}$ and k_{2}) are listed in Table 1.
In contrast to the previous studies, where the free empirical parameters ${P}^{*},e,{k}_{\mathrm{T}}$ and k_{2} 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}$, σ_{3}=σ_{xy}, ${\dot{\mathit{\u03f5}}}_{\mathrm{1}}={\partial}_{x}u+{\partial}_{y}v$, ${\dot{\mathit{\u03f5}}}_{\mathrm{2}}={\partial}_{x}u{\partial}_{y}v$, ${\dot{\mathit{\u03f5}}}_{\mathrm{3}}={\partial}_{y}u+{\partial}_{x}v$, bulk viscosity $\mathit{\zeta}={P}_{\mathrm{p}}(\mathrm{1}+{k}_{\mathrm{T}})/\mathrm{2}{\mathrm{\Delta}}^{\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\ast}$, and replacement pressure $P={P}_{\mathrm{p}}\mathrm{\Delta}/{\mathrm{\Delta}}^{\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\ast}$, Eq. (2) can be split into three decoupled relationships:
The first equation is obtained by taking the trace of Eq. (2). Subtracting the equation for σ_{yy} from the one for σ_{xx} in the set of four relationships (2) yields Eq. (11), while Eq. (12) is just the equation for σ_{xy} extracted from the same set. Equations (10)–(12) are advanced in time s using the Euler scheme with the subcycling time step $\mathit{\delta}{t}_{\mathrm{s}}=\mathit{\epsilon}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{d}}$ (ε=0.00347; see Table 1):
Hereinafter, all the fields without superscripts are taken from the previous time steps. After that, velocity components in Eq. (1) are advanced in time:
Here $m=\stackrel{\mathrm{\u0303}}{\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:
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 Bgrid. 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 twostage 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 socalled strong constraint statespace formulation of the problem, which minimizes a userdefined cost function on the manifold whose structure is specified by the equations of the model. The cost function 𝒥 was defined by
Here, W and $\stackrel{\mathrm{\u0303}}{W}$ denote the nonzero elements of the userdefined (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 gridscale 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
with respect to $\mathit{X},\phantom{\rule{0.125em}{0ex}}\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},\phantom{\rule{0.125em}{0ex}}\mathit{\delta}\widehat{\mathit{X}}$ and δC) to zero:
As it is seen, Eq. (22) simply presents the nonlinear model trajectory specified by a given set of control variables. The second relationship contains the derivatives of 𝒥 and $\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 limitedmemory quasiNewtonian (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:

specify a control vector $\mathit{C}=\mathit{\left\{}\mathit{u}\right(\mathrm{0},\mathit{x}),\phantom{\rule{0.125em}{0ex}}A(\mathrm{0},\mathit{x}),\phantom{\rule{0.125em}{0ex}}h(\mathrm{0},\mathit{x}),\phantom{\rule{0.125em}{0ex}}{\mathit{C}}_{\mathrm{p}}(t,\mathit{x}\left)\mathit{\right\}}$ at the nth iteration;

run the forward model (Eqs. 13–19) to compute X_{n}, the derivatives δ𝒥∕δX and the value of 𝒥_{n} (Eq. 20);

run the adjoint model forced by $\mathit{\delta}\mathcal{J}/\mathit{\delta}\mathit{X}$ (cf. Eq. 23) to obtain ${\widehat{\mathit{X}}}_{n}$;

compute 𝒥_{m}∕δC (Eq. 24);

update the control variables using the M1QN3 software;

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 finitedifference TL and adjoint models are completely defined by the finitedifference 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 finitedifference 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 firstorder 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 firstorder 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 wellknown 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 smallscale 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 ridgelike 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},\phantom{\rule{0.25em}{0ex}}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{\u03f5}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:
, 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.
It is noteworthy that the behavior of Φ(ϵ) in similar experiments with the 1D VP sea ice model and corresponding TLA models (Stroh et al., 2019) agreed well with the ϵ^{2} dependence of the Taylor expansion (blue line in Fig. 1). We speculate that, similar to the 1D case, stability criteria of the 2D VP system may also be weakly affected by the TL transition due to the specific nature of the linearization in the implicit momentum equation solver (see Appendix B). These considerations may indicate an additional attractive feature of the VP rheology in the practical 4DVar applications. Note also that, to the best of our knowledge, the existing 4DVar 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_{2}. 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 CryoSat2, 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 CryoSat2 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 ICESat2 satellite provides highresolution (40 m) freeboard estimates (https://icesat2.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 highresolution SIV with an accuracy of 0.005 m/s (Komarov and Barber, 2014), suggesting a possibility of highprecision 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 highresolution 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 statespace 4DVar 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).
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_{2} and k_{T}. The control vector in these experiments included only parameters k_{2} and k_{T}, and firstguess initial conditions for SIT and SIC were assumed to be errorfree and were not adjusted. Observations of the sea ice velocities, thickness and concentration were assimilated. In the second group of experiments forced by gyreshaped winds (GYRE0, GYREW), 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 firstguess 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.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.
The firstguess solution was forced by the same wind but with k_{T}=0. Figure 2c shows the firstguess 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 4DVar 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.
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 sparsegrid representation of RPs in the 4DVar experiments. If the initial SIT and SIC distributions are not errorfree, 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_{2}
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{\u0303}}{k}}_{\mathrm{1}}$ = 8, the basal stress parameter ${\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}}_{\mathrm{b}}$ = 20 and the maximum basal stress parameter k_{2} = 15 N/m^{3}. 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 firstguess solution had the same initial conditions, wind forcing and RP distribution with the exception that k_{2} 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 firstguess solution.
The variational assimilation of SIV, SIT and SIC observations, targeted at optimization of k_{2}, 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_{2} is shown in Fig. 3g. The maximum values of k_{2} are very close to the true k_{2} = 15 N/m^{3}. Note that the true k_{2} value was specified ad hoc and the grounding effect formally could be reached with smaller values of k_{2}. This is clearly demonstrated in Fig. 3g, where k_{2} is about 12 N/m^{3} 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_{2} was added to the cost function, and the optimized field of k_{2} 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_{2}, which should invoke faster convergence of the minimization procedure.
To demonstrate the robustness of the k_{2} optimization with respect to possible variations in ${\stackrel{\mathrm{\u0303}}{k}}_{\mathrm{1}}$, we conducted a similar experiment with a smaller true ${\stackrel{\mathrm{\u0303}}{k}}_{\mathrm{1}}$ = 2.5. As seen in Fig. 4a and b, the decrease in ${\stackrel{\mathrm{\u0303}}{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 firstguess solution for this experiment was the same as in Fig. 3c with k_{2} = 0, i.e., all ice moving eastward. However, after 4DVar DA and optimization of k_{2}, the reconstructed solution (Fig. 4d) is nearly identical to the true solution (Figure 4b). This optimized k_{2} (Fig. 4c) has the same spatial structure but a slightly smaller (∼ 14 N/m^{2}) maximum value compared to the experiment where ${\stackrel{\mathrm{\u0303}}{k}}_{\mathrm{1}}$ = 8.
The parameterization of grounding landfast ice also includes the critical thickness parameter ${\stackrel{\mathrm{\u0303}}{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{\u0303}}{k}}_{\mathrm{1}}$ than k_{2} (Lemieux et al., 2015) because ${\stackrel{\mathrm{\u0303}}{k}}_{\mathrm{1}}$ can be interpreted as a scaling coefficient in the definition of the critical ice thickness ${\stackrel{\mathrm{\u0303}}{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_{1} optimization through the gradientbased 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_{2}.
4.1 Cyclonic gyre experiments (GYRE0/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^{2} 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 Gaussianshaped 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^{2}. 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^{2}) 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 biasfree 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^{2} for SIC, SIT, SIV and wind stress, respectively. The initial conditions for the firstguess 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^{2}, e = 2 and true wind forcing.
Despite the exact wind forcing, the firstguess 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}},\phantom{\rule{0.25em}{0ex}}{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}},\phantom{\rule{0.25em}{0ex}}{\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_{0} 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
where Ω is the model domain and index k = 1,2 enumerates suboptimal optimization stages: k=1 using the initial conditions control vector C_{0} 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).
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 firstguess 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 firstguess 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^{2}, 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_{0} from the control vector C_{0}. 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 GYRE0/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^{2}. 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).
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^{2}. 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 firstguess initial SIT/SIV conditions and data for the PIZ OSSE were derived from the true solution in a similar way as for the GYRE0 OSSE. Similarly, we specify the first guess with P^{*} = 27.5 kN/m^{2}, e = 2 and the exact wind forcing. The optimization was conducted in three steps: by first optimizing C_{0}, then optimizing C_{rh}, and, finally, simultaneously optimizing C_{0} 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 GYRE0/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.
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^{2}, or about 30 %–35 % of the domainaverage 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 onehalf 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^{2} (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^{2}. 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.
The presented study continues our previous efforts (Stroh et al., 2019) and addresses the feasibility of retrieving spatially varying RPs through the 4DVar 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_{1}), initial conditions and wind forcing for a singlecategory 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 diffusionbased 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 4DVar 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 Jacobianfree Krylov solvers for VP rheology. This approach was employed in the ocean model by Nechaev et al. (2005), where the internal matrixfree biconjugate 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 Jacobianfree 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 4DVar 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_{2} 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_{2} and k_{1} requires a very small (∼ 10) number of iterations, which suggests the possibility of efficiently incorporating their optimization into a panArctic 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_{1}, which confines k_{2} variability to shallow areas and narrow straits areas in a highresolution panArctic sea ice model. It should be noted that retrieval of k_{1} is not straightforward because of nondifferentiability of the grounding parameterization scheme (Eq. 8). Optimizing k_{1} 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_{1}, 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 rootmeansquare 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 CryoSat2 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^{2}) 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 4DVar applications utilized a relatively small assimilation window (3 d) with an eye towards improving shortterm 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 4DVar reconstruction derived from a model driven by incorrect winds. In this case, both the model simulation and the 4DVar results will be inaccurate. However, a properly formulated 4DVar 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_{2} and k_{T} remained virtually unchanged. The incoming satellite platforms (e.g., ICESat2, https://icesat2.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 elastobrittle rheology of Dansereau et al., 2016), may be within the focus of our studies in the near future.
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{\u03f5}}}_{\mathrm{1},\mathrm{2},\mathrm{3}},m,\mathit{\zeta},{P}_{\mathrm{p}}$, P and C_{b} of X are given by
Hereinafter, all the terms with variations of the RPs and atmospheric forcing and variations of the quantities that may contain those are underlined. Taking Eqs. (A1)–(A7) into account, the TL equations of Eqs. (13)–(19) are given by
The gradients with respect to RPs and atmospheric forcing control variables are given by
It is important to note that stability of the adjoint model, being strictly related to the stability of the TL model, contains the Newtonian damping given by the boldfaced terms $\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 4DVar 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).
The VP system of equations is obtained by eliminating the time derivative in Eq. (2), resolving the remainder with respect to σ,
and substituting (B1) into the momentum equation:
where we used the notation P for the tensor in the rhs of Eq. (B1):
The nonlinear system of VP Eqs. (B1)–(B4) is solved in two stages. First, Eq. (B2) is stepped forward using a semiimplicit scheme, which takes the fields of P and Δ^{∗} from the previous time step, treating only $\dot{\mathit{\u03f5}}$ semiimplicitly. 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 4DVar 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.
All inputs to the model experiments are from idealized conditions as specified in the paper. The numerical scheme is provided in large detail in the Appendix and can be easily replicated. Model results in this study are securely stored at the Navy DSRC archive server and can be accessed after obtaining an account at the facility. The corresponding author can be contacted for information to access the archived data once an account has been established.
All authors provided substantial contribution to the models' development, interpretation of the results and writing the manuscript.
The authors declare that they have no conflict of interest.
This research was funded in part under the Naval Research Laboratory’s 6.2 project Modeling Arctic Landfast Ice (program element 62435N) and the 6.1 project Optimization Rheology and Advancing Sea Ice Forecast (program element 61153N). Oceana Francis was supported by the Coastal Hydraulics Engineering Resilience (CHER) Lab, the Civil and Environmental Engineering Department, and the Sea Grant College Program at the University of Hawai`i at Mānoa. Special thanks to the three anonymous reviewers for their comments, which have significantly improved this paper.
This research has been supported by the Naval Research Laboratory’s 6.2 project Modeling Arctic Landfast Ice (program element 62435N) and the 6.1 project Optimization Rheology and Advancing Sea Ice Forecast (program element 61153N) funded by the Office of Naval Research.
This paper was edited by Christian Haas and reviewed by three anonymous referees.
Alexandrov, V., Sandven, S., Wahlin, J., and Johannessen, O. M.: The relation between sea ice thickness and freeboard in the Arctic, The Cryosphere, 4, 373–380, https://doi.org/10.5194/tc43732010, 2010.
Anderson, D. L. and Weeks, W. F.: A theoretical analysis of sea icestrength, EOS T. Am. Geophys. Un., 30, 632–640, 1958.
Arnold Jr., C. P. and Dey, C. H.: Observing System Simulation experiments: Past, Present and Future, B. Am. Meteorol. Soc., 67, 687–695, https://doi.org/10.1175/15200477(1986)067, 1986.
Bouillon, S., Fichefet, T., Legat, V., and Madec, G.: The elastic–viscous–plastic method revisited, Ocean Model., 71, 2–12, 2013.
Cummings, J. A. and Smedstad, O. M.: Variational data assimilation for the global ocean, in: Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications, SpringerVerlag, Berlin Heidelberg, 303–343, https://doi.org/10.1007/9783642350887_13, 2013.
Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elastobrittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359, https://doi.org/10.5194/tc1013392016, 2016.
Errico, R. M.: What is an Adjoint model?, B. Am. Meteorol. Soc., 78, 2577–2591, 1997.
Giering, R. and Kaminski, T.: Recipes for adjoint code construction, ACM T. Math. Software, 24, 437–474, 1998.
Gilbert, J. C. and Lemaréchal, C.: Some numerical experiments with variablestorage quasiNewton algorithms, Math. Progr., 45, 407–435, https://doi.org/10.1007/BF01589113, 1989.
Goldberg, D. E.: Genetic Algorithms in Search, Optimization and Machine Learning, AddisonWesley, Reading, MA, 432 pp., 1989.
Goldberg, D. N. and Heimbach, P.: Parameter and state estimation with a timedependent adjoint marine ice sheet model, The Cryosphere, 7, 1659–1678, https://doi.org/10.5194/tc716592013, 2013.
Gray, J. M. N. T. and Killworth, P. D.: Stability of viscousplastic sea ice rheology, J. Phys. Oceanogr., 25, 971–978, 1995.
Harder, M. and Fischer, H.: Sea ice dynamics in the weddell sea simulated with an optimized model, J. Geophys. Res.Oceans, 104, 11151–11162, 1999.
Heimbach, P., Menemenlis, D., Losch, M., Campin, J. M., and Hill, C.: On the formulation of seaice models. Part 2: Lessons from multiyear adjoint sea ice export sensitivities through the Canadian Arctic Archipelago, Ocean Model., 33, 145–158, https://doi.org/10.1016/j.ocemod.2010.02.002, 2010.
Hibler, W.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979.
Hibler, W. D. and Walsh, J. E.: On modeling seasonal and interannual fluctuations of arctic sea ice, J. Phys. Oceanogr., 12, 1514–1523, 1982.
Hoteit, I., Cornuelle, B., Kohl, A., and Stammer, D.: Treating strong adjoint sensitivities in tropical eddypermitting variational data assimilation, Q. J. Roy. Meteor. Soc., 131, 3659–3682, 2005.
Hunke, E. and Dukowicz, J.: An elasticviscousplasticmodel for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997.
Hunke, E. C.: Viscousplastic sea ice dynamics with the EVP model: Linearization issues, J. Comput. Phys., 170, 18–38, 2001.
Hunke, E. C. and Lipscomb, W. H.: CICE: The Los Alamos sea ice model. Documentation and software user's manual version 4.0, Tech. Rep. LACC06012, Los Alamos Natl. Lab., Los Alamos, NM, 2008.
Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos Sea Ice Model documentation and software users manual version 4.1, Tech. Rep. LACC06012, T3 Fluid Dynamics Group, Los Alamos Natl. Lab., Los Alamos, NM, 2010.
Ingber, L.: Very fast simulated reannealing, Math. Comput. Model., 12, 967–973, 1989.
Juricke, S., Lemke, P., Timmermann, R., and Rackow, T.: Effects of stochastic ice strength perturbation on Arctic finite element sea ice modeling, J. Climate, 26, 3785–3802, https://doi.org/10.1175/JCLID1200388.1, 2013.
Kauker, F., Kaminski, T., Karcher, M., Giering, R. , Gerdes, R. and Voßbeck, M.: Adjoint analysis of the 2007 all time Arctic seaice minimum, Geophys. Res. Lett., 36, L03707, https://doi.org/10.1029/2008GL036323, 2009.
Kimmritz, M., Danilov, S., and Losch, M.: The adaptive EVP method for solving the sea ice momentum equation, Ocean Model., 101, 59–67, https://doi.org/10.1016/j.ocemod.2016.03.004, 2016.
Koldunov, N. V., Aizinger, V., Rakowsky, N., Scholz, P., Sidorenko, D., Danilov, S., and Jung, T.: Scalability and some optimization of the FinitevolumE Sea ice–Ocean Model, Version 2.0 (FESOM2), Geosci. Model Dev., 12, 3991–4012, https://doi.org/10.5194/gmd1239912019, 2019.
Komarov, A. S. and Barber, D. G.: Sea ice motion tracking from sequential dualpolarization RADARSAT2 images, IEEE T. Geos. Remote S., 52, 121–136, 2014.
König Beatty, C. and Holland, D. M.: Modeling landfast sea ice by adding tensile strength, J. Phys. Oceanogr., 40, 185–198, https://doi.org/10.1175/2009JPO4105.1, 2010.
Kreyscher, M., Harder, M., and Lemke, P.: First results of the SeaIce Model Intercomparison Project (SIMIP), Ann. Glaciol., 25, 8–11, 1997.
Kreyscher, M., Harder, M., Lemke,P., and Flato, G. M.: Results of the sea ice model intercomparison project: Evaluation of sea ice rheology schemes for use in climate simulations, J. Geophys. Res.Oceans, 105, 11299–11320, 2000.
Laxon, S. W., Giles, K. A., Ridout, A. L., Wingham, D. J., Willatt, R., Cullen, R., Kwok, R., Schweiger, A., Zhang, J., Haas, C., Hendricks, S., Krishfield, R., Kurtz, N., Farrell, S., and Davidson, M.: CryoSat‐2 estimates of Arctic sea ice thickness and volume Geophys. Res Lett., 40, 732–737, https://doi.org/10.1002/grl.50193, 2013.
Le Dimet, F. X.: A general formalism of variational analysis, CIMMS report 22, Cooperative Institute for Mesoscale Meteorological Studies, Norman, OK, USA, 34 pp., 1982.
Le Dimet, F. X. and Talagrand, O.: Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects, Tellus, 38A, 97–110, 1986.
Lemieux, J.F. and Tremblay, B.: Numerical convergence of viscous‐plastic sea ice models, J. Geophys. Res.Oceans, 114, C05009, https://doi.org/10.1029/2008JC005017, 2009.
Lemieux, J.F., Tremblay, B., Thomas, S., Sedlacek, J., and Mysak, L. A.: Using the preconditioned generalized minimum residual (GMRES) method to solve the seaice momentum equation, J. Geophys. Res.Oceans, 113, C10004, https://doi.org/10.1029/2007JC004680, 2008.
Lemieux, J.F., Knoll, D., Tremblay, B., Holland, D., and Losch, M.: A comparison of the Jacobianfree NewtonKrylov method and the EVP model for solving the sea ice momentum equation with a viscousplastic formulation: a serial algorithm study, J. Comp. Phys., 231, 5926–5944, 2012.
Lemieux, J.F., Knoll, D., Losch, M., and Girard, C.: A secondorder accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics, J. Comput. Phys., 263, 375–392, 2014.
Lemieux, J.F., Tremblay, L. B., Dupont, F., Plante, M., Smith, G. C., and Dumont, D.: A basal stress parameterization for modeling landfast ice, J. Geophys. Res.Oceans, 120, 3157–3173, https://doi.org/10.1002/2014JC010678, 2015.
Lemieux, J.F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.Oceans, 121 , 7354–7368, 2016.
Lewis, J. M. and Derber, J. C.: The use of adjoint equations to solve a variational adjustment problem with advective constraints, Tellus A, 37, 309–322, https://doi.org/10.3402/tellusa.v37i4.11675, 1985.
Losch, M., Fuchs, A., Lemieux, J. F., and Vanselow, A.: A parallel Jacobianfree NewtonKrylov solver for a coupled sea ice ocean model, J. Comp. Phys., 257, 901–911, 2014.
Massonnet, F., Goosse, H., Fichefet, T., and Counillon, F.: Calibration of sea ice dynamic parameters in an oceansea ice model using an ensemble Kalman filter, J. Geophys. Res.Oceans, 119, 4168–4184, https://doi.org/10.1002/2013JC009705, 2014.
Massonnet, F., Fichefet, T., and Goosse, H.: Prospects for improved seasonal Arctic sea ice predictions from multivariate data assimilation, Ocean Model., 88, 16–25, 2015.
Metzger, E. J., Smedstad, O. M., Thoppil, P. G., Hurlburt, H. E., Cummings, J. A., Wallcraft, A. J., Zamudio, L., Franklin, D. S., Posey, P. G., Phelps, M. W., Hogan, P. J., Bub, F. L., and DeHaan, C. J.: US Navy operational global ocean and Arctic ice prediction systems, Oceanography, 27,32–43, https://doi.org/10.5670/oceanog.2014.66, 2014.
Miller, P. A., Laxon, S. W., Feltham, D. L., and Cresswell, D. J.: Optimization of a sea ice model using basinwide observations of Arctic sea ice thickness, extent, and velocity, J. Climate, 19, 1090–1108, 2006.
Nechaev, D. A., Panteleev, G., and Yaremchuk, M.: Reconstruction of the Circulation In Limited Regions of an Ocean With Open Boundaries: Climatic Circulation In the Tsushima Strait, Oceanology, 45, 761–780, 2005.
Nguyen, A. T., Menemenlis, D., and Kwok, R.: Arctic iceocean simulation with optimized model parameters: Approach and assessment, J. Geophys. Res., 116, C04025, https://doi.org/10.1029/2010JC006573, 2011.
Nichols, N. K.: Data Assimilation: Aims and Basic Concepts, in: Data Assimilation for the Earth System, NATO Science Series (Series IV: Earth and Environmental Sciences), edited by: Swinbank, R., Shutyaev, V., and Lahoz, W. A., Springer, Dordrecht, 2003.
Nichols, N. K.: Mathematical concepts of data assimilation, in: Data assimilation: making sense of observations, edited by: Lahoz, W., Khattatov, B., and Menard, R., Springer, Dordrecht, 2010.
Nicolsky, D. J, Romanovsky, V. E., and Panteleev, G. G.: Estimation of soil thermal properties using insitu temperature measurements in the active layer and permafrost, Cold Reg. Sci. Technol., 55 120–129, 2009.
Nitta, T.: Some analyses of observing systems simulation experiments in relation to the First GARP Global Experiment, GARP Working Group on Numerical Experimentation, Report No 10, 1–35, Plan for U.S. Participation in the Global Atmospheric Research Program, National Academy of Sciences, Washington, DC, USA, 1969.
Panteleev, G., Rogers, W. E., Yaremchuk, M., Shen, H., Rainville, L., and Grout, J.: Floe Size Mapping from Satellite SAR Images and Icewatch Observations in the Beaufort Sea during Autumn 2015, Tech. Rep. NRL/MR/7322199903, Naval Researh Laboratory, Stennis Space Center, Mississippi, USA, 2019.
Penenko, V. V.: Methods of Numerical Simulation of Atmospheric processes, Gidrometeoizdat, Lenigrad, 350 pp., 1981.
Posey, P. G., Metzger, E. J., Wallcraft, A. J., Preller, R. H., Smedstad, O. M., and Phelps, M. W.: Validation of the 1/128 Arctic Cap Nowcast/Forecast System (ACNFS), Tech. Rep. NRL/MR/7320109287, Naval Res. Lab., Stennis Space Center, Mississippi, USA, 2010.
Ryan, P. A. and Münchow, A.: Sea ice draft observations in Nares Strait from 2003 to 2012, J. Geophys. Res., 122, 3057–3080, https://doi.org/10.1002/2016JC011966, 2017.
Schwegmann, S., Haas, C., Fowler, C., and Gerdes, R.: A comparison of satellitederived seaice motion with driftingbuoy data in the Weddell Sea, Antarctica, Ann. Glaciol., 52, 103–110, 2011.
Simmonds, I. and Rudeva, I.: The great Arctic cyclone of August 2012, Geophys. Res. Lett., 39, L23709, https://doi.org/10.1029/2012GL054259, 2012.
Stroh, J. N, Panteleev, G., Yaremchuk, M., Francis, O., and Allard, R.: Toward optimization of rheology in sea ice models through data assimilation, J. Atm. Oceanic Tech., 36, 2365–2382, https://doi.org/10.1175/JTECHD180239.1, 2019.
Sumata, H., Kwok, R., Gerdes, R., Kauker, F., and Karcher, M.: Uncertainty of arctic summer ice drift assessed by highresolution SAR data, J. Geophys. Res.Oceans, 120, 5285–5301, 2015.
Sumata, H., Kauker, F., Karcher, M., and Gerdes, R.: Simultaneous Parameter Optimization of an Arctic Sea Ice–Ocean Model by a Genetic Algorithm, Mon. Weather Rev., 147, 1899–1926, https://doi.org/10.1175/MWRD180360.1, 2019.
Thorndike, A. S. and Colony, R.: Sea ice motion in response to geostrophic winds, J. Geophys. Res., 87, 5845–5852, https://doi.org/10.1029/JC087iC08p05845, 1982.
Tilling, R., Ridout, A., and Shepher, A.: Estimating Arctic sea ice thickness and volume using CryoSat2 radar altimeter data, Adv. Space Res., 62, 1203–1225, 2018.
Toyota, T. and Kimura, N.: An examination of the sea ice rheology for seasonal ice zones based on ice drift and thickness observations, J. Geophys. Res.Oceans, 123, 1406–1428, 2018.
Tremblay, L. B. and Mysak, L. A.: Modeling sea ice as a granular material, including the dilatancy effect, J. Phys. Oceanogr., 27, 2342–2360, 1997.
Tremblay, L. B., and Hakakian, M.: Estimating the sea ice compressive strength from satellite derived sea ice drift and NCEP reanalysis data, J. Phys. Oceanogr., 36, 2165–2172, 2006.
Tschudi, M., Meier, W., Stewart, J., Fowler, C., and Maslanik, J.: Polar Pathfinder daily 25 km EASEGrid sea ice motion vectors, version 4, dataset 0116, NASA NationalSnow and Ice Data Center Distributed Active Archive Center, Boulder, CO, USA, https://doi.org/10.5067/INAWUWO7QH7B, 2019.
Uotila, P., Farrell, S. O., Marsland, S., and Bi, D.: A seaice sensitivity study with a global oceanice model, Ocean Model., 51, 1–18, https://doi.org/10.1016/j.ocemod.2012.04.002., 2012.
Vancoppenolle, M., Fichefet, T., Goosse, H., Bouillon, S., Madec, G., and Maqueda, M. A. M.: Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 1. Model description and validation, Ocean Model., 27, 33–53, 2009.
Wunsch, C.: The Ocean Circulation Inverse Problem, Cambridge Univ. Press, Cambridge, 442 pp., 1996.
Yaremchuk, M. and Martin, P.: On Sensitivity Analysis within the 4DVAR Framework, Mon. Weather Rev., 142, 774–787, 2014.
Yaremchuk, M. and Nechaev, D.: Covariance localization with diffusionbased correlation models, Mon. Weather Rev., 141, 848–860, 2013.
Yaremchuk, M., Nechaev, D., and Panteleev, G.: A method of successive corrections of the control subspace in the reducedorder variational data assimilation, Mon. Weather Rev., 137, 2966–2978, 2009.
Yaremchuk, M., Townsend, T., Panteleev, G., Hebert, D., and Allard, R.: Advancing short‐term forecasts of ice conditions in the Beaufort Sea, J. Geophys. Res., 124, 807–820, https://doi.org/10.1029/2018JC014581, 2019.
Zhang J. and Hibler III, W. D.: On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res., 102, 8691–8702, 1997.
Zhang, J. and Rothrock, D.: Modeling global sea ice with a thickness and enthalpy distribution model in generalized curvilinear coordinates, Mon. Weather Rev., 131, 845–861, 2003.
Zhang, Y.F. and Bitz, C. M.: Insights on sea ice data assimilation from perfect model observing system simulation experiments, J. Climate, 31, 5911–5926, 2018.
 Abstract
 Introduction
 Sea ice model and its 4DVar implementation
 Optimization of the landfast parameters
 Optimization of the ice strength and axes ratio fields
 Summary and discussion
 Appendix A: Tangent linear numerics
 Appendix B: On the stability of TL models with VP rheology
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Sea ice model and its 4DVar implementation
 Optimization of the landfast parameters
 Optimization of the ice strength and axes ratio fields
 Summary and discussion
 Appendix A: Tangent linear numerics
 Appendix B: On the stability of TL models with VP rheology
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References