Parameter Optimization in Sea Ice Models with Elastic-Viscoplastic Rheology

. Ice rheology formulation is a key component of modern sea ice modeling. 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 ﬁxed parameters empirically adjusted to optimize the model performance. In this study, we consider a spatially variable representation of these parameters to extend a two-dimensional EVP sea ice model with a formulation similar to CICE6. Feasibility of 5 optimization of the rheological and landfast sea ice parameters is assessed via idealized variational data assimilation experiments with synthetic observations of ice concentration, thickness and velocity. Taking into account that sea ice observations are available daily, the experiments are conﬁgured for a 3-day data assimilation window in a rectangular basin. It is found that the tangent linear and adjoint models featuring EVP rheology are unstable, but can be stabilized by adding a Newtonian damping term into the adjoint equations. The set of the observation system simulation experiments shows that landfast parameter 10 distributions can be reconstructed after 5-10 iterations of the minimization procedure. Optimization of the sea ice initial conditions and spatially varying parameters in the equation for the stress tensor requires more computation, but provides a better hindcast of the sea ice state and the internal stress tensor. Analysis of the inaccuracy in the wind forcing and errors in the sea ice thickness observations have shown reasonable robustness of the variational DA approach and feasibility of its application to available and incoming observations.


Introduction
Due to significant decline of the sea ice volume in the Arctic Ocean (AO) in the last decades and increase of the maritime activity in the AO, the role of an accurate sea ice hindcast/forecast became an important component of the regional data assimilation systems. Currently, there are several community sea ice models broadly used for modeling and reconstruction of 20 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 visco-plastic (VP) rheology proposed by Hibler (1979). In the last decades, several numerical approaches have been proposed for solving the VP problem (Hibler, 1979;Zhang and Hibler, 1997;Lemieux et al, 2008). These approaches are based on implicit solvers and require a significant number of iterations to achieve full convergence. Application of more efficient quasi-Newtonian solvers suffers 25 from the lack of robustness and these are usually applied in the sea ice models of intermediate resolution (Lemieux et al., 2012;Losch et al., 2014). In addition, implicit VP solvers are less convenient for implementing on massively parallel supercomputer architectures. Despite these inconveniences the currently known four-dimensional variational (4dVar) sea ice models utilize an implicit VP solver for time integration of the tangent linear and adjoint (TLA) models (Kauker et al, 2009;Heimbach et al, 2010). 30 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 widely used 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 Con-35 sortium. It is widely used in sea ice and coupled sea-ice modeling (e.g. Posey et al, 2010, Metzger et al, 2014. There are multiple examples of 3dVar, ensemble-based DA systems utilizing the CICE (e.g., Zhang and Bitz, 2018), but to our knowledge there is no 4dVar DA system based on the CICE model yet.
In the majority of sea ice models, the sea ice rheology is controlled by three major parameters (P ⇤ , e, and ↵) describing respectively the dimensional maximum ice strength per unit thickness, the ratio of yield ellipse major axes, and non-dimensional 40 scaling of ice strength with its compactness. For modeling landfast ice near the coast and in the straits (so-called arching phenomena), an additional parameter k T has been introduced to model the tensile strength (Konig 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 , ↵ b , for better parameterization of the landfast ice grounded in the shallow regions. Typically, the above mentioned rheological parameters (RPs) are constants and their values are 45 defined empirically from multiple numerical experiments. RPs such as P ⇤ and e reflect the model parameterization rather than physics and thus are directly unobservable (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 the sea ice drift were diagnosed to vary within 27.5 kN/m 2 (Hibler III and Walsh, 1982), 15-20 kN/m 2 (Kreyscher et al., 1997(Kreyscher et al., , 2000, and 30-45 kN/m 2 (Tremblay and Hakakian, 2006). These studies indicate the existence of significant variations of the P ⇤ estimates, which may be attributed 50 to both non-physical considerations (such as spatially variable model resolution), and spatio-temporal variations of the Arctic sea ice.
Numerical experiments of Lemieux et al., (2016) with a coarse resolution pan-Arctic CICE-NEMO model have shown that k T = 0.2 provides the best agreement with landfast strength observations in the Kara Sea, when the ellipse axes ratio ranges within 1.2-1.4. The most sensitive parameters for sea ice grounding are the critical thickness parameter k 1 and maximum basal 55 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 the landfast modeling for different parts of the Arctic Ocean, suggesting that these parameters are the functions of local environmental conditions. Physical properties of the sea ice also strongly depend on sea ice salinity, temperature and ice age (e.g. Anderson and Weeks, 1958) indicating that different categories of sea ice may have different rheological properties. 60 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 "trialand-error" approach, involving multiple runs of a sea ice model with different sets of RPs (e.g. Miller et al, 2006, Uotila et al, 2012, while others utilized the more advanced methods based on the Green's function approach (Nguyen et al, 2011), ensemble Kalman filtering (Massonnet et et al., 2013) and genetic algorithms (Sumata et al., 2019). However, all of these attempts sought a relatively small set of spatially invariant sea ice model parameters which provided 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 the need of extremely high computational resources.
Also note, that the above listed algorithms are not suitable for simultaneous optimization of other model paratemers such as 70 initial and open boundary conditions, or external forcing fields.
The major objective of our study is to find a numerically feasible method for simultaneous optimization of multiple model parameters including the spatially varying RPs, initial conditions and forcing fields in the sea ice models based on EVP solvers (e.g. CICE model). We suggest that, if successful, this approach can be adopted for optimizing the RPs in the operational sea ice models (e.g. Posey et al., (2010), Metzger et al, (2014) and provide more accurate short term (3-7 days) sea ice forecast.

75
Feasibility of reconstructing the spatially varying fields of P ⇤ and e (as well as other model parameters) through the variational DA of synthetic observations of Sea Ice Velocity/Concentration/Thickness (SIV/SIC/SIT) was recently analyzed by Stroh et al, (2019) in the framework of the 1d (zonal) VP sea ice model. It was found that variational DA allows reasonable reconstruction of spatially varying P ⇤ and e in the regions with strong convergence and significantly improves short range hindcast/forecast of the sea ice state. In particular, it was shown that optimization of spatially varying P ⇤ and e provides more 80 accurate reconstructions of the ridging areas, that cannot be achieved by optimizing initial sea ice state only.
Note, that optimization of the RPs through the 4Dvar DA approach allows us to efficiently use all available sea ice observations including sea ice velocity, that are rarely used for assimilation in sea ice DA systems. The latter is due to weak sensitivity of the sea ice state with respect of the ice velocity (e.g. Kauker, et al., 2009). Roughly speaking, the 4dVar DA approach allows us to use sea ice velocity observations for adjustment of the RPs and/or atmospheric forcing in an appropriate manner resulting 85 in 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 2d sea ice model based on the EVP rheology formulation of Lemieux et al., (2016). Our approach is based on application of the 4Dvar and Observing System Simulation Experiment (OSSE) concepts (e.g. Nichols 2003Nichols , 2010 and follow the conventional twin-data experiment approach (e.g. Goldberg and Heimbach, 2013). Similarly to Stroh et al, (2019), we developed the corresponding 90 EVP TLA models and analyzed the feasibility of the optimization of spatially varying ellipse ratio and sea ice compressive strength. In addition, we also analyzed the feasibility of optimizing two of the landfast sea ice parameters introduced by Lemieux et al.(2016). Through the 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 optimization of the basal stress parameters in the ice sheets model (Goldberg and Heimbach, 2013). 95 Currently, satellite sea ice observations are typically available daily with a reasonably dense spatial resolution. Analysis of the SAR images (e.g. Panteleev et al., 2019) indicates that in the marginal sea ice zone, the pancake/cake ice with floe sizes of ⇠1-20 m may be easily replaced by floes exceeding 1 km in size in one week. As a consequence, we configured the OSSEs with a 3-day DA window assuming that such approach should have more impact on short term sea ice forecast.
The paper is organized as follows: Section 2 describes the implemented sea ice model, the details of the TLA codes and 100 generation of synthetic observations and the first guess solution used in OSSE experiments. Results of these experiments are described in Sections 3 (optimization landfast sea ice parameters) and 4 (optimization of compressive strength and ellipse ratio), with 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.

105
This section provides details of the sea ice model formulation, its associated linearizations, outlines the variational data assimilation system used for optimization of the 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.

110
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 SI 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, tr are the divergence and trace operators, k is the vertical unit vector, I is the 2⇥ 2 identity matrix, h, A, and u = {u, v} are the 2d fields of ice thickness, concentration and velocity, and and✏ are the 2d fields of ice stress and the deformation rate tensors: The scalar field used for normalization of the rhs in (2) is computed using the following expression (e.g., Hunke, 2001): To avoid numerical singularities at✏ = 0, the values of are limited from below by the additional parameter˜ ⇤ = 10 10 s 1 , so that ⇤ (✏) = max( ,˜ ⇤ ).

125
The empirical parametersT d , P ⇤ , k T and e in equations (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 Hibler's (1979) rheology: The typical values of the ice strength parameter P ⇤ and↵ are listed in Table 1.

130
The bottom and ocean stresses in eq. (1) were parameterized in accordance with Lemieux et al. (2015; 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 135 angle ⇥ counterclockwise,ũ 0 = 10 5 m/s, and u w is the water velocity (set to zero in the present study). The values of other parameters (C w , ⇢ w ,↵ b ,k 1 and k 2 ) are listed in Table 1.
In contrast to the previous studies, where the free empirical parameters P ⇤ , e, k 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 C ⌘ {h, A, u} using variational DA technique. (2) with respect to ice velocities and RPs. Note that the first term in the rhs of the linearized eq. (2) is proportional to the first derivatives of the velocity perturbations u. As a consequence, the components of are linear in the first derivatives of u after taking the explicit time step t s in the linearized eq. (2). Moreover, the first-order derivatives in u keep their presence in the rhs of the linearized eq. (1) due to spatial variability of the background fields in equations (6) and (7).

155
This property of the TL equations of the subsystem (1)-(2) may require additional care in the specification of 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 non-linear 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 non-linear model. In particular, the non-linear stability criterion could be violated, for example, in the areas of strong ice 160 convergence where the background field is characterized by the large values of SIT gradients which may cause, as it follows from (7), large coefficients before the first-order derivatives of u in the TL code of 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.
observed similar instability of the TL EVP solver in the MIT sea ice model (personal communication).
Instability of the linearized codes in the strongly non-linear regimes of the parent model is a well known phenomenon in the ocean general circulation models (OGCMs). A heuristic solution of the problem was proposed by Hoteit et al, (2004) who added extra diffusion in TLA codes to suppress unstable small-scale harmonics. This kind of treatment is achieved, however, at the expense of the reduction of the TLA code accuracy (e.g., Yaremchuk et al, (2009)). Later, Yaremchuk and Martin,170 (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 with 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 175 ridge-like structures (i.e. having a wide spatial spectrum of the SIT component) in the areas of ice convergence. As a consequence they cannot be efficiently damped using an isotropic diffusion tensor in the TL momentum equations. A straightforward solution is to introduce spatially inhomogeneous diffusion tensor field (e.g. Yaremchuk and Nechaev, 2013), with local anisotropy derived from the background solution. However, this would require a considerable reduction of t s due to very large diffusion along the ridges. As a simple alternative, we used a more simple Newtonian friction term, which homogeneously 180 damps the entire spectrum of small perturbations. Numerical experimentation has shown that this approach worked generally well using Newtonian damping time scaleT d of 0.7 t. Additional experiments have shown thatT d has to be decreased to 3 t d in the cases of stronger sea ice convergence in the regions with thick (h > 3 m) ice.
Testing the validity of the stabilized TLA codes was done in a way similar to Yaremchuk et al, (2009). The initial conditions for the model thickness and concentration fields C = {h(x, 0), A(x, 0)} were slightly perturbed by the realizations of a random 185 function R (viz. C(x, 0) ! C(x, 0) + ✏R(x) ⌘ C + C), and the model and its TL version were integrated for t = 5 t. After that, the dependence of the normalized difference between the non-linear solution and its TL approximation was checked by computing the following quantity: where N is the state vector of the non-linear model at time t generated by the initial conditions listed in the argument, and | · | 190 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 behavior of (✏) in similar experiments with 1d VP sea ice model and corresponding TLA models 195 (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 In the OSSE experiments described below, the control variables included 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., (2015Lemieux et al., ( , 2016 and listed in Table 1.

Simulated observations and cost functions 205
In all OSSE experiments 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 pre-processing, 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 provides high resolution (40m) freeboard estimates https://icesat-2.gsfc.nasa.gov/ 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  and temporal decorrelation scale of 7 days. Taking into account that high resolution satellite SIC, SIV are currently available on a daily basis over the entire Arctic Ocean and assuming they can be interpolated within every daily time frame, we set synthetic observations to be available in all the space-time grid points (600⇥8990 ⇡ 5.4 · 10 6 ) of the model domain.
The standard state-space 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  Table 1).

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 of the 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 first guess initial conditions for SIT and SIC were assumed to be error-free and were not adjusted. Observations of the sea ice velocities, 240 thickness and concentration were assimilated. In the second group of experiments forced by gyre-shaped winds (GYRE-0, GYRE-W), we analyzed the feasibility of optimizing P ⇤ and e in the regions with spatially and temporally varying SIT and SIC. In the third group, we analyzed the feasibility of optimizing P ⇤ and e in the Pack Ice Zone (PIZ) where SIT varies in space and SIC is equal to 1. We also explore the impact of the P ⇤ and e optimization 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 245 the first guess sea ice initial conditions were accordingly disturbed. A list OSSEs and with short descriptions are assembled in

Formation of landfast ice in the deep narrow straits and between islands is a well known phenomenon in the Canadian
Archipelagos and in the Kara Sea (e.g. Lemieux et al, 2016). In the Nares Strait, landfast ice is observed periodically and 255 typically its boundary has an arching shape (e.g. Ryan and Munchow, 2017).
To mimic this phenomenon, the sea ice model was configured in a narrow zonal domain forced by steady zonal wind for of k T was set to 0.6. Figure 2b shows, that after 3 days 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 265 increasing SIC along the southern boundary.
The first guess solution was forced by the same wind but with k T =0. Figure 2c shows the first guess state at the end of the assimilation window (t=3 days). 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 10km wide polynya developing at the western boundary. The landfast ice is completely absent in this solution due to the absence of tensile strength in 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 Figure 2d. The optimized k T is shown in Figure 2f and 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 the 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 275 function, which is proportional to the integral of k 2 T 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 (Figure 2e), which is the result of the small dimension of the k T control vector due to representation of the RPs on a sparse grid in 4dVar experiments. If the initial SIT and SIC distributions are not error-free, it is also possible to optimize the initial SIT and SIC in addition to k T . This kind of the optimization provides better SIT/SIC hindcast but requires more iterations to find the optimal 280 solution (dashed line in Figure 2e).

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 (Figure 3) with zonally varying depth 285 ranging between 3m at the western boundary and 33 m at the eastern boundary (Figure 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.25m and 2.5m as shown in Figure 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 parameterk 1 =8, basal stress parameter↵ b =20, and the maximum basal stress parameter k 2 =15 290 N/m 3 . Figure   The variational assimilation of the SIV, SIT and SIC observations targeted at optimization of k 2 demonstrated a significant improvement of SIT (not shown) and SIC, clearly seen in Figure 3d. In particular, the optimized solution produced the landfast ice region and the polynya location that are nearly identical to those observed in true solution (Figure 3b).
The optimized field of k 2 is shown in Figure 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 305 is clearly demonstrated in Figure 3g, where k 2 is about 12 N/m 3 over the major part of the landfast ice area in true solution (cf. 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) (Figure 3f). Note also, that according to equations (1) and (8), sea ice acceleration is directly proportional to k 2 which should invoke faster convergence of the 310 minimization procedure.
To demonstrate the robustness of the k 2 optimization with respect to possible variations ink 1 we conducted a similar experiment with smaller truek 1 =2.5. As it seen from Figure 4a,b, the decrease ofk 1 causes a considerable decrease of the landfast area in the southwestern corner. The landfast ice polynya also moved about 100km closer to the western coastline and is now confined to the shallower region. The first guess solution for this experiment was the same as in Figure 3c with k 2 =0, 315 i.e. all ice moved eastward. However, after 4dVar DA and optimization of k 2 the reconstructed solution (Figure 4d) is nearly identical to the true solution (Figure 4b). The optimized k 2 is shown at Figure 4c and has the same spatial structure but slightly smaller (⇠ 14 N/m 2 ) maximum value as compared to the experiment withk 1 =8.
The parametrization of the grounding landfast ice also includes the critical thickness parameterk 1 , which was kept fixed in the described experiments. According to multiple numerical simulations, the total landfast ice area is more sensitive to 320 variations ofk 1 than k 2 (Lemieux et al, 2015) becausek 1 can be interpreted as a scaling coefficient in the definition of the critical ice thicknessk 1 h c = Ah b (cf. eq. 8), and works as a "switch" which defines the areas of potential landfast ice generation.
However, the discontinuity of the Heaviside step function present in equation 8 significantly complicates k 1 optimization through the gradient-based variational method. Formally, this problem can be regularized (e.g., Nicolsky et al., 2009), but this approach requires additional optimization of the additional regularization parameter. In practice, this approach could be 325 computationally less efficient. In the light of this consideration, we limit our study with analysis of the feasibility of the landfast ice parameterizations to optimizing k T and k 2 .  Figure 5e,f within the following ranges: 22.5 kN/m 2  P ⇤  32.5 kN/m 2 and 1  e 3. These ranges were adopted from various studies (e.g. Hibler and Walsh, 1982;Kreyscher et al. 1997Kreyscher et al. , 2000Tremblay and Hakakian 2006;Lemieux et al. 2016). The true wind forcing had a form of the Gaussian-shaped cyclone with stationary position whose 335 strength gradually increased by 1.5 times during the 3-day assimilation window. The resulting wind stress at t=0 is shown in Figure 5c and had a maximum value of 0.7 N/m 2 . Initial SIV conditions were determined by a 100 minute 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 days under the applied atmospheric forcing. Figures 5c-d show the trace of the stress tensor P tr = tr /2 at the beginning of the true state integration and after 3 days. The P tr distribution has a clear maximum near the 340 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.5m) 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.5m) sea ice.
Noisy SIC, SIT and SIV observations were generated by adding spatially and temporally correlated noise (with the corre-345 lation scales of 150 km and 7 days) 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 section 2.3, i.e. they have similar absolute errors as most of the real observations available, nowadays. In the experiments we did not introduce any bias to ice observations since it is a common assumption in the existing DA systems.
The magnitudes of the imposed noise correspond to the errors of the respective observational data sets, with the amplitudes   ducted an additional optimization of the full control vector C = {C 0 , C rh }. 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 non-linear optimization problems whose cost functions may have multiple minima. Figure 6a-d compares 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 sub-optimal optimization stages: k = 1 using the initial conditions control vector C 0 only, and k = 2 employing of the full control vector C = {C 0 , C rh }. It was found that S u reduced almost 1.5 times after the additional optimization of the rheological parameters C rh (Figure 6a,b). Visual comparison of the sub-optimal and fully optimized SIC shows a certain improvement after C rh as well. For example, the local minimum of the sub-optimal SIC in the region with coordinates 700-1500 km and 180 km is about 0.7 (white contours in Figure 6a), while the fully 365 optimized SIC has a minimum of 0.6 which agrees perfectly with the true SIC distribution (white contours in Figure 5b).
In this experiment we found that optimization of C rh yields only a marginal improvement of the SIT distribution. In particular, std(h k opt h 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 substantial difference between the first guess and observed SITs. Another possible reason is that C rh is not well controlled by the initial 370 SIT distribution at the time scale of the relatively short 3-day assimilation window.
As expected, optimization of the rheological parameters C rh provides the major impact on the correction to the internal stress tensor. Figures 6c,d, show that standard deviation of the differences between the true and the optimized values of P tr decreased by 35% after the additional optimization of the RPs. Note also, that the maximum value of fully optimized P tr ( Figure 6d) is located in the close proximity of the true P tr maximum (Figure 5d), while the maximum in sub-optimal P tr is 375 shifted almost 400 km northward. Comparison with true P ⇤ demonstrates that the optimized field agrees well with true P ⇤ almost everywhere with the only exception observed in the southwest corner of the domain (Figure 6e). This is caused by the substantial SIV divergence (Figure 6a) which diminishes the role of rheological forcing in the region.
The reconstruction of e is not as accurate, as P ⇤ with the extreme values of the optimized e being 1.2 and 2.8, while the respective true values are 1 and 3. Note, however, that spatial locations of extrema agrees well with true distribution. The 380 exception is the local minimum in the south-western corner, where both e and P ⇤ disagree with true values being close to their first guess values. Note, that significant improvement of P tr and the internal stress components discussed above are directly related to the more accurate optimization of the RP and SIV fields because both P ⇤ and e control the structure of the stress tensor in eq. (2).
To analyze impact of the inaccuracy in the wind forcing we conducted an additional experiment where the center of the 385 cyclonic disturbance was displaced 90 km westward mimicking a systematic error in the hypothetical atmospheric forecast. i.e. only by 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.5kN, i.e. about 40-50% worse than in the experiment with exact forcing, but it is important to note that the maximum in the P tr is still in very good agreement with the true solution.

395
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 maxima in the western and eastern parts of the region and a minimum in the center of the domain (Figure 7c), while the minimum of the e in the western part demonstrates certain agreement with true e distribution (Figure 5e,f). Note, that inaccurate wind forcing affects the accuracy of P ⇤ retrievals to a less degree than e.

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 days of integration (Figure 5a,b). Due to the exponential dependence of P on 1 A, (eq. 7) internal stress decreases exp(4) ⇡ 50 times, and therefore has a minor rheological impact of 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 405 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 .
The model domain, initial SIT, and P ⇤ and e distributions were the same as in GYRE-0/W-OSSEs. However, unlike the cyclonic wind of the previous experiment, the atmospheric forcing is applied as a 20 m/s eastward wind at the western boundary 410 which reverses zonal direction across the breadth of the domain (Figure 8c). In time, the wind speed was amplified 1.5 times (to 30 m/s) by the end of the DA window. The resulting wind stress at t=3d is shown in Figure 7c and has a maximum amplitude about 0.5 N/m 2 .
The temporal evolution of the true SIT and SIV is shown at Figure 8a The first guess initial SIT/SIV conditions and data for the PIZ-OSSE were derived from the true solution in a similar 420 way as for GYRE-0-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 C rh , and, finally, simultaneous optimization of C 0 and C rh . Figure 9a,c shows SIT, and the difference between optimized and true SIV at t=3d after optimization of the initial conditions only. Interestingly, despite less spatial variations in wind forcing, optimization of the initial conditions does not allow for 425 accurate reconstruction of the SIV as we observed in the GYRE-0/W-OSSEs. The maximum errors in the eastern part of the domain are about 0.1 m/s, being comparable in the 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 sub-optimal distribution of P tr (Figure 9c) differs significantly from true P tr (Figure 8d) both quantitatively and qualitatively. For example std(P 1 tr (opt)-P tr (true)) is 11.6 kN/m 2 , or about 30-35% of the domain-average value of P tr . The qual- itative difference is probably more important because sub-optimal P tr distribution fails to provide two maximums discussed above. Instead, the sub-optimal distribution of P tr has only one maximum located in the center (Figure 8d).
Additional optimization of the rheological parameters C rh significantly improved the reconstructed SIV practically everywhere, so, the formal measure of the uncertainty S u decreased almost two times from 56 m/s down to 30 m/s (Figure 9b).
Similar improvements are visible in P tr distribution (Figure 9d). The std[P 2 tr (opt) P tr (true)] decreased to 7.4 kN/m 2 (by 435 40%) and the fully optimized P tr has two maxima as in the true P tr distribution (Figure 8d).
The optimized P ⇤ and e are shown at Figure 9e,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. The effect could probably be attributed to the region with zero convergence along the western boundary where the rheology does not play a significant role. There is also some quantitative difference between optimized and true e in 440 the central part of the model domain but, qualitatively, the reconstructed field of e has all features of the true distribution. The presented study continues our previous efforts (Stroh et al, 2019) and addresses the feasibility of retrieving spatiallyvarying RPs through the 4dVar assimilation the satellite observations of SIV, SIC and SIT. To do the analysis, we developed TLA codes with respect of all rheological parameters (except k 1 ), initial conditions, and wind forcing for a single-category sea 445 ice model recently proposed by Lemieux et al., (2016). The dynamical core of this model is based on conventional formulation of the EVP rheology (Hunke and Dukowicz, 1997;Hunke and Lipscomb, 2008) and parameterizations of the grounding and arching land fast ice recently proposed by Lemieux et al., (2015Lemieux et al., ( , 2016 and Koning 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 450 employed the Lax-Wendroff scheme for ice advection (Appendix A).
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 455 ones observed in the case of diffusion-based stabilization, and thus the Newtonian scheme can be successfully used in sea ice models based on the EVP solvers. In the last decade several modifications of the EVP solvers were proposed to improve the convergence (e.g., Bouillon et al, 2013, Kimmritz et al, 2016, Koldunov et al, 2019. We assume that these ideas could benefit 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 more efficient stabilization of the EVP TLA 460 codes. On the other hand, simple analysis ( Figure 1, dashed line) indicates that the TLA model for the 1D VP sea ice solver (Stroh et al, 2019) is stable provided that 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), who 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 465 (e.g. Hibler et al. (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 Jacobian-free Krylov solvers for VP rheoplogy. This approach was employed in the ocean model by Nechaev et al., (2005), where the internal matrix-Free BiCG solver was constructed using the model's TLA codes and in the IMEX algorithm, where the preconditioned Flexible GMRES algorithm was implemented 470 in the Jacobian-free mode (Lemieux 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 reasonable reconstruction of the RPs. The numerical experiments included two series of the 4dVar DA runs.
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 the landfast ice simulations in various parts of the Arctic Ocean (Lemieux et al., 2015), especially in the shallow sea and the narrow straits of the Canadian Archipelago. In this study, these parameters were specified by spatially variable functions with a reduced number (10-36, Section 2.4) of free parameters that were optimized to fit surface observations. The conducted OSSEs demonstrate that spatially varying landfast ice parameters k T and k 2 , responsible 480 for grounding and arching phenomena, can be optimized at a relatively low computational cost (5-12 iterations).
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 high-resolution pan-Arctic sea ice model. It should be noted that retrieval of k 1 is not straightforward because of non-differentiability of the grounding parameterization scheme (eq. 8). Optimizing k 1 in sea ice 485 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 which is about two times smaller than the errors of the CryoSat-2 observations. The OSSEs also assessed sensitivity of the results with respect to systematic errors in wind forcing.

495
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 1d sea ice model featuring VP parameterization of the RPs, supporting a higher impact of the spatially varying P ⇤ than e on the DA quality. We also observed that regions of less accurate P ⇤ reconstructions are typically collocated in the regions of strong sea ice divergence and/or SIC concentrations below 0.8, i.e. with the regions where 500 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 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 very important because it better informs avoidance of regions with excessive compressive stress.
The OSSE with zonal wind convergence in the pack ice (A=1) demonstrated even better quality in reconstructing e and 505 especially P ⇤ . This can be attributed to a stronger impact of the 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 .
We utilized a relatively small assimilation window (3 days) keeping in mind 4dVar application for improving short-term sea ice forecasts in the ice pack and ice edge zones. In these regions, ice rheology typically changes at the time scales of several days 510 (e.g. Panteleev et al, 2019) due to variations in the dynamic and/or thermodynamic forcing. In particular, wind variability may cause profound changes of the ice edge, while in the ice pack, wind forcing is the major driving factor of polynya formation, due to ice divergence and the ridging processes, when rheological forcing becomes important. Finally, sea ice observations are typically available on a regular (daily) basis, in time, 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 515 to the complex structure and uncertainties of the respective observation operators. These developments in data acquisition and preprocessing will result in improved observability of the RP fields within the data assimilation systems.
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 520 k 2 and k T remained virtually unchanged. RPs estimation on the local scale through the 4dVar DA approach presented here. It is also noteworthy, that the extension of the similar 4dVar DA based on the VP rheology may improve the accuracy of the reconstructed rheological parameters due to better stability of the TLA codes in the VP formulation. This extension, as well as more realistic inversions employing more complex rheological hypotheses (e.g., the Maxwell elasto-brittle rheology of Dansereau et al., (2016)), may be within the focus of our studies in the near future.

530
Appendix A: Model numerics and the cost function In deriving the TL and adjoint models, we closely followed the EVP numerics given in Lemieux et al. (2016).
Introducing notations 1 = xx + yy , 2 = xx yy , 3 = xy , ✏ 1 = @ x u + @ y v, ✏ 2 = @ x u @ y v, ✏ 3 = @ y u + @ x v, bulk viscosity ⇣ = P p (1 + k T )/2 ⇤ , and replacement pressure P = P p / ⇤ , equation (2) can be split into three decoupled relationships: The first equation (A1) is obtained by taking the trace of (2). Subtracting the equation for yy from the one for xx in the set of four relationships (2) yields (A2), while (A3) is just the equation for xy extracted from the same set. Equations (A1-A3) 540 are advanced in time s using the Euler scheme with the subcycling time step t s = "T d (" = 0.00347, see Table 1): Hereinafter all the fields without superscripts are taken from the previous time steps. After that, velocity components in equation

545
(1) are advanced in time: Here m =⇢hA and ⌧ a is the atmospheric forcing. After the equations (A4-A8) are advanced in time for 400 time steps, the ice thickness and concentration fields are updated using the obtained velocity and the Lax-Wendorff scheme.
All spatial derivatives present in (A4-A10) were discretized by finite differences on the Arakawa B-grid. At the boundaries we used the condition u = 0. Initial conditions for u were set either to zero or defined through the model integration for 1 hour.
Initial conditions for A and h were specified as arbitrary functions.

555
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 MIT 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. TAMS (Giering and Kaminski, 1998;560 http://autodiff.com/tamc/), OpenAD (https://www.mcs.anl.gov/OpenAD/), we briefly outline the details of the development of the EVP TLA models.
Perturbations of the auxiliary functions ✏ 1,2,3 , m, ⇣, P p , P, 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 (A11-A17) into account, the TL equations of (A4-A10) are given by The last two terms in (A23-A24) are TL for Lax-Wendroff non-linear diffusion, andD is the Laplacian operator. In the 585 variational assimilation experiments we used the strong constraint state-space formulation of the problem, minimizing the cost function: Here, W ,W denotes the non-zero elements of the (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 590 ⌦.
The optimization problem under the constraints imposed by the sea ice model (eq. A4-A10) requires minimization of the modified cost function whereX t are the fields of Lagrangian multipliers containing the adjoints of the model state variables X t specified in all 595 the grid points in ⌦. The sequence X t depends on control fields C = [C 0 , C p ] which include the initial state of the model C 0 ⌘ X 0 , and other control fields C p which contain rheological parameters, and atmospheric forcing errors. The latter two were defined on a sparser grid (in every 15th node of the computational grid for most of the conducted experiments) and spatially interpolated on the model grid using the bilinear interpolation operator I. The minimization of the cost function with with the additional range constraints for the selected control fields (Section 2.3) performed after each iteration. The M1QN3 algorithm requires a procedure for computation of the cost function gradient J / C with respect to C for a given value of the control fields. This procedure was implemented in two steps.
First, the sequence of model states X t (C) was computed by specifying the first guess values of the control fields C and integrating equations (A4-A10) over the assimilation window for 3 days. In the course of integration, the derivatives J / X t 605 of the cost function (A25) where also computed and stored in memory. Note that this integration can be represented as setting to zero the derivatives of J m with respect to the adjoint fieldsX t in eq. (A26).
On the second sweep, the gradients J m / C are calculated. In the process, J m / X t are set to zero for all t except t = 0 by integrating the adjoint of the TL equations backward in time with the rhs forcing terms equal to J / X t computed on the first (forward) sweep. The adjoint equations, being linear inX t , are obtained by transposition of the TL system matrix in the 610 right hand side of (A18-A24). During the backward integration, all the underlined terms in (A18-A24) containing variations of C p were set to zero to obtain the gradient of J with respect to X 0 ( J m / C 0 ) at the end of integration. In paralell, the gradients J / C p were calculated by computing the respective derivatives of the second termJ in (A26). Numerically, this was done by replacing C p in (A18-A24) by the components ofX (1 + k T )Ah exp( ↵(1 A))(✏ 1˜ 1 + ✏ 2˜ 2 + ✏ 3˜ 3 /2)/ ⇤ (A27) J m e = "I T X t ⇢ 2e (1 + "e 2 ) 2 [ 2˜ 2 + 3˜ 3 + "P p (1 + k T )(✏ 2˜ 2 + ✏ 3˜ 3 /2)/ ⇤ ] P p (1 k T ) e 3 (1 + ") (✏ 2 2 + ✏ 2 3 )˜ 1 + + P p e 3 ⇤2 ✓( ˜ ⇤ )(✏ 2 2 + ✏ 2 3 )


(1 + k T ) ✓ ✏ 1˜ 1 1 + " + ✏ 2˜ 2 + ✏ 2˜ 3 /2 1 + e 2 " ◆ + (1 k T ) It is important to note that stability of the adjont model, being strictly related to the stability of the TL model, contains the Newtonian damping given by the terms " 1,2,3 in eq. (A18-A20). The adjoint model is integrated backward in time and requires the solution of the forward model for the each time step of the backward sub-cycling procedure. This can be 625 achieved either through re-calculating 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. 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 system matrix in (A18-A24).

630
Appendix B: On the stability of TL models with VP rheology The VP system of equations is obtained by eliminating the time derivative in eq. (2), resolving the remainder with respect to = 1 e 2  P I 1 e 2 2 trP , and substituting (B1) into the momentum equation: where we used the notation P for the tensor in the rhs of (B1): The non-linear system of VP equations (B1-B4) is solved in two stages. First, equation (B2) is stepped forward using a semi-implicit scheme, which takes the fields of P and ⇤ from the previous time step, treating only✏ semi-implicitly. On the second stage, equations (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 procedure would not provide the exact derivative of the non-linear scheme with respect 645 to u because equation (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, 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 performace of the above mentioned "incomplete linearization" of the VP model in the 650 4dVar applications, as soon as the stability criterion of the non-linear model (B1-B4) is satisfied. We performed an additional experiment with 1d VP forward and TL models using the modified procedure featuring ten applications of the GMRES solver and ten 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.