Toward a method for downscaling sea ice pressure for navigation purposes

Sea ice pressure poses great risk for navigation; it can lead to ship besetting and damages. Contemporary large-scale sea ice forecasting systems can predict the evolution of sea ice pressure. There is, however, a mismatch between the spatial resolution of these systems (a few kilometres) and the typical dimensions of ships (a few tens of metres) navigating in ice-covered regions. In this paper, the downscaling of sea ice pressure from the kilometre-scale to scales relevant for ships is investigated by conducting high-resolution idealized numerical experiments with a viscous-plastic sea ice model. Results show that sub-grid-scale pressure values can be significantly larger than the large-scale pressure (up to ∼ 4 times larger in our numerical experiments). High pressure at the sub-grid scale is associated with the presence of defects (e.g. a lead). Numerical experiments show significant stress concentration on both sides of a ship beset in sea ice, especially at the back. The magnitude of the stress concentration increases with the length of the lead (or channel) behind the ship and decreases as sea ice consolidates by either thermodynamical growth or mechanical closing. These results also highlight the difficulty of forecasting, for navigation applications, the small-scale distribution of pressure, and especially the largest values as the important parameters (i.e. the length of the lead behind the ship and the thickness of the refrozen ice) are not well constrained.

1 Introduction

With the growing shipping activities in the Arctic and surrounding seas, there is a need for user-relevant sea ice forecasts and products at multiple timescales and spatial scales. An important forecast field for navigation is the internal sea ice pressure (simply referred to as pressure for the rest of this paper). In compact ice conditions, high-pressure events can complicate navigation activities and even pose great risk for ship besetting.

Figure 1 A 24 h forecast of the sea ice pressure (kN m⁻¹) and of the surface winds (m s⁻¹) from the Canadian Arctic Prediction System (CAPS). The forecast was initiated at 00:00 UTC on 29 April 2020. Almost all of the domain is shown in (a), while (b) is a subset of the domain located in the region of Svalbard (the sub-region is defined by the blue rectangle in a). Note that the colour scale is not the same for the two panels.

By solving equations for the momentum balance and for the ice thickness distribution, sea ice models are able to predict the evolution of the pressure field. However, even for high-resolution operational forecasting systems with spatial resolutions of a few kilometres (e.g. Dupont et al., 2015; Hebert et al., 2015), there is a clear mismatch in the spatial scales considered. Indeed, the forecast pressure from the model, which represents the average pressure for a grid cell of a few kilometres wide, is not necessarily relevant for a much smaller ship; there are larger pressure values than the average pressure provided by the sea ice forecasting system (Kubat et al., 2010, 2012; Leisti et al., 2011). Figure 1 shows an example of a pressure forecast from a large-scale forecasting system. The Canadian Arctic Prediction System (CAPS) is a fully coupled atmosphere–sea ice–ocean system developed and maintained by the Canadian Centre for Meteorological and Environmental Prediction. Its domain covers the Arctic Ocean, the North Atlantic and the North Pacific. The spatial resolution of the atmospheric model is ∼ 3 km, while the spatial resolution in this region for the sea ice and ocean models is ∼ 4.5 km. Looking at a specific region that is north of Svalbard (Fig. 1b), it can be observed that the surface winds push the ice toward the coast and create large pressure.

Some researchers have done case studies of compressive or besetting events using large-scale sea ice forecasting systems (e.g. Kubat et al., 2012, 2013; Leisti et al., 2011). These besetting events are all associated with heavy ice conditions. The investigations of Kubat et al. show the importance of the coast on pressure conditions; the sea ice pressure often increases toward the coast.

Mussells et al. (2017) used ship logs and satellite imagery to relate besetting events and density of sea ice ridges. They indeed found that the ship was often beset in areas and times of the year with high ridge densities. Probabilistic models for ship performance in sea ice and likelihood of besetting events have also been developed (e.g. Montewka et al., 2015; Turnbull et al., 2019). Turnbull et al. (2019) argue that the primary cause of the besetting events they studied were the relatively large ice floes encountered by the vessel.

There is also a vast literature on the performance of ships navigating in ice-infested waters and on the estimation of ice resistance, that is the longitudinal forces applied on the ship by the ice (e.g. Lindqvist, 1989; Su et al., 2010; Jeong et al., 2017). These calculations are important for ship design and for operational considerations. Lindqvist (1989) introduced a simple empirical formulation to calculate ice resistance based on a ship's characteristics and ice physical parameters. When sea ice pressure is not considered, the resistance encountered by a ship only depends on processes such as crushing, breaking and displacement of ice floes by the ship's hull. Based on laboratory experiments in ice tanks, Kulaots et al. (2013) extended this empirical approach by also considering the effect of compression on the performance of ships navigating in ice-infested waters. There are also some numerical studies of ice loads on ships transiting in ice-infested waters where sea ice is represented using discrete elements (i.e. the floes; Metrikin and Løset, 2013; Daley et al., 2014) or as a continuum (e.g. Kubat et al., 2010).

In contrast with studies mentioned in the last paragraph, we focus on ship besetting rather than on a ship progressing in an ice-covered region. We also study the downscaling of sea ice pressure from the kilometre scale to scales relevant for navigation activities (tens of metres). Note that this was briefly investigated by Kubat et al. (2010) for a ship transiting through a loose sea ice cover. Kubat et al. (2010) showed that the pressure on the hull of the ship can be 2 orders of magnitude larger than the large-scale pressure. For our numerical experiments, we use a continuum-based viscous-plastic sea ice model. In a first set of simulations, we study how the small-scale pressure depends on the stresses applied at the boundaries, on the ice conditions and on the rheology parameters. The second part of the results is dedicated to shipping applications; we investigate the small-scale pressure field in the vicinity of an idealized ship beset in heavy ice conditions and under compressive stresses. Idealized sea ice modelling studies with a continuum-based approach have been conducted by specifying strain rates at the boundaries (e.g. Kubat et al., 2010; Ringeisen et al., 2019) or by specifying wind patterns (e.g. Hutchings et al., 2005; Heorton et al., 2018). However, to our knowledge, it is the first time that internal stresses are specified at the boundaries.

This paper is structured as follows. In Sect. 2, the sea ice momentum equation and the viscous-plastic rheology are described. The sea ice model used for the numerical experiments is presented in Sect. 3. The approach for prescribing sea ice stresses at the boundaries is presented in Sect. 4. The validation of our experimental setup is done in Sect. 5. The main results are given in Sect. 6. Finally, concluding remarks are provided in Sect. 7.

2 Sea ice momentum equation and rheology

The large-scale sea ice forecasting system solves the sea ice momentum given by

$$\begin{array}{}\text{(1)}& \mathit{\rho }h{\displaystyle \frac{D\mathit{u}}{Dt}}=-\mathit{\rho }hf\widehat{\mathit{z}}\times \mathit{u}+{\mathit{\tau }}_{\mathrm{a}}-{\mathit{\tau }}_{\mathrm{w}}+\mathbf{\nabla }\cdot \mathit{\sigma }-\mathit{\rho }hg\mathbf{\nabla }{H}_{\mathrm{d}},\end{array}$$

where ρ is the density of the ice; h is the ice volume per unit area (sometimes referred to as the mean thickness); $\frac{D}{Dt}$ is the total derivative; f is the Coriolis parameter; $\mathit{u}=u\widehat{\mathit{x}}+v\widehat{\mathit{y}}$ is the horizontal sea ice velocity vector; $\widehat{\mathit{x}}$, $\widehat{\mathit{y}}$ and $\widehat{\mathit{z}}$ are unit vectors aligned with the x, y and z axis of our Cartesian coordinates; τ_a is the wind stress, τ_w is the water stress; σ is the internal ice stress tensor with components σ_ij acting in the jth direction on a plane perpendicular to the ith direction; g is the gravitational acceleration; and H_d is the sea surface height. This two-dimensional formulation, which is obtained by integrating along the vertical, is justified when the ratio between the horizontal and vertical scales of the problem is large (i.e. a ratio of at least 10:1; Coon et al., 1974).

The sea ice pressure is by definition the average of the normal stresses, that is

$$\begin{array}{}\text{(2)}& p=-({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}})/\mathrm{2},\end{array}$$

with a negative sign because, by convention, stresses in compression are negative. The pressure is the first stress invariant (i.e. it does not vary with the choice of the coordinate system). The second stress invariant (q), that is the maximum shear stress at a point, is defined by

$$\begin{array}{}\text{(3)}& q=\sqrt{{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}+{\left[{\displaystyle \frac{({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}})}{\mathrm{2}}}\right]}^{\mathrm{2}}}.\end{array}$$

As the stresses are vertically integrated, the stresses and stress invariants are 2D fields with units of newtons per metre. Because the sea ice stresses are written as a function of the sea ice velocity (see details below), one also obtains the pressure p and the maximum shear stress q when solving the momentum equation for u. Hence, by solving the momentum equation for the large-scale sea ice model, the pressure at every grid point is obtained (we refer to this pressure field as the large-scale pressure).

Here, we consider a small area of sea ice (the size of a grid cell) to be under compressive stresses. The idea is to apply the large-scale pressure at the boundaries of this small area and to simulate the sub-grid-scale sea ice pressure (referred to as the small-scale pressure). We assume here that the ice is neither moving nor deforming (e.g. it is being held against a coast). To further simplify the problem, the wind stress, the water stress, the advection of momentum and the sea surface tilt term are neglected. We wish to find, inside this small domain, the steady-state solution of $\mathit{\rho }h\partial \mathit{u}/\partial t=\mathbf{\nabla }\cdot \mathit{\sigma }$, which is equivalent to finding the solution of $\mathbf{\nabla }\cdot \mathit{\sigma }=\mathrm{0}$. The stresses are modelled according to the viscous-plastic rheology with an elliptical yield curve (Hibler, 1979). With this rheology, the relation between the stresses and the strain rates can be written as

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{ϵ}}}_{ij}+[\mathit{\zeta }-\mathit{\eta }]{\dot{\mathit{ϵ}}}_{kk}{\mathit{\delta }}_{ij}-P{\mathit{\delta }}_{ij}/\mathrm{2},i,j=\mathrm{1},\mathrm{2},\end{array}$$

where δ_ij is the Kronecker delta; ${\dot{\mathit{ϵ}}}_{ij}$ are the strain rates defined by ${\dot{\mathit{ϵ}}}_{\mathrm{11}}=\frac{\partial u}{\partial x}$, ${\dot{\mathit{ϵ}}}_{\mathrm{22}}=\frac{\partial v}{\partial y}$ and ${\dot{\mathit{ϵ}}}_{\mathrm{12}}=\frac{\mathrm{1}}{\mathrm{2}}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})$, ${\dot{\mathit{ϵ}}}_{kk}={\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}$; ζ is the bulk viscosity; η is the shear viscosity; and P is a term which is a function of the ice strength.

The bulk and shear viscosities are, respectively,

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{P_{\mathrm{p}}}{\mathrm{2}\mathrm{△}}},\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{\eta }=\mathit{\zeta }{e}^{-\mathrm{2}},\end{array}$$

where P_p is the ice strength;
$\mathrm{△}={\left[({\dot{\mathit{ϵ}}}_{\mathrm{11}}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}^{\mathrm{2}})(\mathrm{1}+e^{-\mathrm{2}})+\mathrm{4}{e}^{-\mathrm{2}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\mathrm{2}}+\mathrm{2}{\dot{\mathit{ϵ}}}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{22}}(\mathrm{1}-e^{-\mathrm{2}})\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}$; and e is the aspect ratio of the ellipse, i.e. the ratio of the long and short axes of the elliptical yield curve.

Following Hibler (1979), the ice strength P_p is parameterized as

$$\begin{array}{}\text{(7)}& P_{\mathrm{p}}=P^{*}h\exp [-C(\mathrm{1}-A\left)\right],\end{array}$$

where P^* is the ice strength parameter; A is the sea ice concentration; and C is the ice concentration parameter, an empirical constant set to 20 (Hibler, 1979) such that the ice strength decreases quickly with the ice concentration. Unless otherwise stated, the rheology parameters P^* and e are, respectively, set to 27.5 kN m⁻² and 2.

When △ tends toward zero, Eqs. (5) and (6) become singular. To avoid this problem, ζ is capped following the approach of Hibler (1979). It is expressed as

$$\begin{array}{}\text{(8)}& \mathit{\zeta }={\displaystyle \frac{P_{\mathrm{p}}}{\mathrm{2}{\mathrm{\Delta }}^{*}}},\end{array}$$

where ${\mathrm{\Delta }}^{*}=max(\mathrm{\Delta },{\mathrm{\Delta }}_{\text{min}})$ with ${\mathrm{\Delta }}_{\text{min}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ s⁻¹.

We use a replacement method similar to the one presented in Kreyscher et al. (2000). The P term in Eq. (4) is given by

$$\begin{array}{}\text{(9)}& P=\mathrm{2}\mathit{\zeta }\mathrm{△}.\end{array}$$

The replacement method is commonly used in sea ice models to prevent unrealistic deformations of the sea ice cover when there is no external forcing.

3 Experimental setup

The McGill sea ice model is used for the numerical experiments. We use revision 333 with some modifications, described below, for specifying stresses at the boundaries.

Considering a domain of a few kilometres by a few kilometres wide (representing a grid cell of a large-scale sea ice forecasting system), the idea is to use the model at very high resolution for studying the distribution of pressure inside that domain. To do so, the model was modified so that internal stresses can be specified at the boundaries (instead of the usual Dirichlet condition, i.e. u=0, at land–ocean boundaries and the Neumann condition at open boundaries with gradients of velocity equal to 0). These stresses at the boundaries represent the integrated effect of the wind and water stresses (like one would get from a large-scale model). The next section gives more details about the implementation of the stress boundary conditions.

For the experiments, the domain is a square of dimensions 5.12 km×5.12 km. It is subdivided into small squared grid cells of dimensions Δx×Δx, with Δx taking one of the following values depending on the experiment: 10, 20, 40, 80, 160, 320, 640 or 1280 m. The size of the domain was chosen because it is close to the average size of CAPS sea ice grid cells and because 5120 m divided by the Δx listed above gives an integer number of small grid cells. For simplicity, we refer to this domain as our 5 km × 5 km domain.

The momentum equation $\mathit{\rho }h\partial \mathit{u}/\partial t=\mathbf{\nabla }\cdot \mathit{\sigma }$ is solved at time levels Δt, 2Δt, 3Δt, …, where Δt is the time step. We introduce the index $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$, which refers to these time levels. Using a backward Euler approach for the time discretization, the momentum equation is written as

$$\begin{array}{}\text{(10)}& \mathit{\rho }h{\displaystyle \frac{({\mathit{u}}^n-{\mathit{u}}^{n-\mathrm{1}})}{\mathrm{\Delta }t}}=\mathbf{\nabla }\cdot {\mathit{\sigma }}^n.\end{array}$$

The spatial discretization of Eq. (10) on the McGill model Arakawa C-grid leads to a system of non-linear equations that is solved using a Jacobian-free Newton–Krylov (JFNK) solver (the most recent version is described in Lemieux et al., 2014). The ice starts from rest. The time step is 30 min. At each time level, the non-linear convergence criterion is reached when the Euclidean norm of the residual has been reduced by a factor of 10. The maximum number of non-linear iterations is set to 500. The steady-state solution is assumed when the maximum velocity difference between two time levels is lower than 10⁻⁹ m s⁻¹. As the ice is assumed to be held against the coast, the simulated velocities are very small (i.e. most of the ice cover is in the viscous regime). Our numerical simulations therefore provide 2D static fields of the internal stresses inside this small domain. Thermodynamic processes and advection of h and A are neglected for all the numerical experiments described in this paper.

For some of the numerical experiments, a digitized ship is placed inside the domain. The digitized ship is simply defined as a rigid body by using land cells. The boundary conditions on the contour of the ship are therefore no outflow and no slip (i.e. u=0), which is the usual Dirichlet approach for land cells. This allows us to investigate the distribution of small-scale pressure around this idealized ship.

4 Boundary conditions for the small domain

The boundary conditions are imposed the same way on the four sides of the small domain. Hence, to shorten the paper, only the treatment on the west side of the domain is explained here. The McGill model uses an Arakawa C-grid; the center of the cell is the point for tracers (e.g. h and A), while the velocity components are positioned on the left side (for u) and lower side (for v). To avoid confusion with the indices i and j for the stresses σ_ij and the strain rates ${\dot{\mathit{ϵ}}}_{ij}$, the indices l and m are, respectively, used to identify the grid cells along the x and y axes. The cell at the southwesternmost location of the domain has indices l=1 and m=1. Figure 2 shows one of the grid cells in the first column of the domain (on the west side). The left side of the grid cell is on the west boundary of the domain. The sides of the domain are referred to as west (W), east (E), south (S) and north (N).

Figure 2 One grid cell on the western boundary of the domain with indices l=1 and m. This figure shows the location of the velocity components on the C-grid of the McGill model. The variables h and A are positioned at the tracer point t. Some variables (e.g. σ₁₂) are also calculated at the node (n) point. The stresses (σ^W and τ^W) applied at the western boundary are shown with purple arrows.

On the west side of the domain, a normal stress (σ^W) and a shear stress (τ^W) are applied. The momentum balance for the u component is comprised of the terms $\partial {\mathit{\sigma }}_{\mathrm{11}}/\partial x$ and $\partial {\mathit{\sigma }}_{\mathrm{12}}/\partial y$. Inside the domain, these terms are approximated by second-order centered differences. At the boundaries, however, a one-sided first-order approximation is employed for $\partial {\mathit{\sigma }}_{\mathrm{11}}/\partial x$. Hence, $\partial {\mathit{\sigma }}_{\mathrm{11}}/\partial x$ at the u location $u(l,m)=u_{\left(lm\right)}$ with l=1 is approximated as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{11}}}{\partial x}}\sim {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}\left(\mathrm{1}m\right)}-{\mathit{\sigma }}_{\left(m\right)}^{\mathrm{W}}}{\mathrm{\Delta }x/\mathrm{2}}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{11}\left(\mathrm{1}m\right)}=[{\mathit{\zeta }}_{\left(\mathrm{1}m\right)}+{\mathit{\eta }}_{\left(\mathrm{1}m\right)}][u_{\left(\mathrm{2}m\right)}-u_{\left(\mathrm{1}m\right)}]\mathrm{\Delta }{x}^{-\mathrm{1}}+[{\mathit{\zeta }}_{\left(\mathrm{1}m\right)}-{\mathit{\eta }}_{\left(\mathrm{1}m\right)}][v_{(\mathrm{1}m+\mathrm{1})}-v_{\left(\mathrm{1}m\right)}]\mathrm{\Delta }{y}^{-\mathrm{1}}-P_{\left(\mathrm{1}m\right)}/\mathrm{2}$ is evaluated at the t point.

On the other hand, the term $\partial {\mathit{\sigma }}_{\mathrm{12}}/\partial y$ only depends on the boundary conditions, that is

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{12}}}{\partial y}}\sim {\displaystyle \frac{{\mathit{\tau }}_{(m+\mathrm{1})}^{\mathrm{W}}-{\mathit{\tau }}_{\left(m\right)}^{\mathrm{W}}}{\mathrm{\Delta }y}}.\end{array}$$

For the v component v_(1m) (which is at a distance of Δx∕2 from the boundary), there is no special treatment for $\partial {\mathit{\sigma }}_{\mathrm{22}}/\partial y$. However, the second-order treatment of the term $\partial {\mathit{\sigma }}_{\mathrm{12}}/\partial x$ follows

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{12}}}{\partial x}}\sim {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{12}\left(\mathrm{2}m\right)}-{\mathit{\tau }}_{\left(m\right)}^{\mathrm{W}}}{\mathrm{\Delta }x}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{12}\left(\mathrm{2}m\right)}={\mathit{\eta }}_{\left(\mathrm{2}m\right)}[u_{\left(\mathrm{2}m\right)}-u_{(\mathrm{2}m-\mathrm{1})}]\mathrm{\Delta }{y}^{-\mathrm{1}}+{\mathit{\eta }}_{\left(\mathrm{2}m\right)}[v_{\left(\mathrm{2}m\right)}-v_{\left(\mathrm{1}m\right)}]\mathrm{\Delta }{x}^{-\mathrm{1}}$ is evaluated at the n point.

Even though u_(1m) is located at the boundary, it is solved along with v_(1m) and all the other velocity components in the domain by the non-linear solver.

In our simulations, ${\mathit{\sigma }}_{\left(m\right)}^{\mathrm{W}}={\mathit{\sigma }}^{\mathrm{W}}$, and ${\mathit{\tau }}_{\left(m\right)}^{\mathrm{W}}={\mathit{\tau }}^{\mathrm{W}}$; i.e. they do not vary with m along the boundary (same idea for the other sides of the domain). Furthermore, for numerical stability (see Appendix A), the normal stress on the east side (σ^E) has to be equal to σ^W. Similarly, σ^S=σ^N and ${\mathit{\tau }}^{\mathrm{W}}={\mathit{\tau }}^{\mathrm{E}}={\mathit{\tau }}^{\mathrm{S}}={\mathit{\tau }}^{\mathrm{N}}$.

5 Model validation

The McGill model has, over the years, been extensively tested (e.g. Lemieux et al., 2014; Bouchat and Tremblay, 2017; Williams and Tremblay, 2018). A few simple experiments were conducted in order to validate the implementation of the new stress boundary conditions.

In all the experiments, normal and shear stresses are applied at the four boundaries of the 5 km × 5 km domain. For a given set of sea ice conditions, the steady-state solution of Eq. (10) is obtained. This provides us with the velocity field defined on the Arakawa C-grid. As the stresses and invariants are functions of the sea ice conditions and velocity (see Eqs. 2–9), static 2D fields of the internal stresses and invariants are easily obtained.

Compared to realistic pan-Arctic simulations, the simplicity of the problem allows one to obtain analytical solutions for specific cases. In a first validating experiment, the thickness (h) and concentration (A) fields are, respectively, set to 2 m and 1 everywhere on the domain. By specifying ${\mathit{\sigma }}^{\mathrm{W}}={\mathit{\sigma }}^{\mathrm{E}}={\mathit{\sigma }}^{\mathrm{S}}={\mathit{\sigma }}^{\mathrm{N}}$ = −10 kN m⁻¹ (i.e. p = 10 kN m⁻¹) and ${\mathit{\tau }}^{\mathrm{W}}={\mathit{\tau }}^{\mathrm{E}}={\mathit{\tau }}^{\mathrm{S}}={\mathit{\tau }}^{\mathrm{N}}$ = 0 kN m⁻¹ at the boundaries, the shear stress should be 0 everywhere inside the domain while the pressure field should be constant and equal to 10 kN m⁻¹. This is indeed what is obtained from the numerical experiment (not shown). With p = 10 kN m⁻¹, a 2 m ice cover is able to resist this compressive stress; that is the ice should be in the viscous regime. Using the definition of the stresses from Eq. (4), we obtain $p=P/\mathrm{2}-\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{I}}$, where ${\dot{\mathit{ϵ}}}_{\mathrm{I}}={\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}$ is the divergence. It is easy to demonstrate that the analytical solution for the divergence is ${\dot{\mathit{ϵ}}}_{\mathrm{I}}=-p{\mathrm{\Delta }}_{\text{min}}/P_{\mathrm{p}}=-\mathrm{3.63}\times {\mathrm{10}}^{-\mathrm{10}}$ s⁻¹. This is exactly what is obtained with the model (not shown).

Figure 3 Pressure field for Δx = 10 m (a) and Δx = 20 m (b). The thickness field is 2 m everywhere except a 1 km long, 40 m wide horizontal lead in the middle of the domain. The normal stresses at the boundaries are −10 kN m⁻¹. Panel (c) shows PDFs of the pressure in the 5 km × 5 km domain for Δx = 10 m (cyan) and Δx = 20 m (dashed magenta). Bins of 0.25 kN m⁻¹ were used to build the PDFs.

We also verify that we obtain the same results when a lead is present within the physical domain for different spatial resolution (Δx). For example, Fig. 3 shows the pressure field for a 1 km long, 40 m wide lead resolved with a Δx of 10 m (Fig. 3a) and for the same lead resolved with Δx = 20 m (Fig. 3b). The thickness of the level ice (h_l) around the lead is 2 m. The maximum pressure at 10 m resolution is 35.79 kN m⁻¹, while the maximum pressure at 20 m is 33.15 kN m⁻¹. From these simulated 2D pressure fields, probability density functions (PDFs) are calculated using bins of 0.25 kN m⁻¹. They are shown in Fig. 3c, which further demonstrates that the simulated fields are very similar at 10 and 20 m resolutions.

The effect of the same lead but oriented differently in the domain was also tested. The PDF of the pressure field is exactly the same whether the lead is oriented horizontally (west–east) or vertically (south–north; not shown). The spatial distribution of pressure is qualitatively the same when orienting the lead diagonally. The PDF of pressure for this diagonal lead is similar to the PDF of the vertical and horizontal ones, although we find that the maximum pressure is usually a bit smaller (not shown). This is likely a consequence of the spatial discretization of a finite width lead on a Cartesian grid.

In a last set of experiments for the validation, we also checked that the presence of relatively nearby boundaries does not affect our conclusions. In the first experiment, with Δx = 20 m, a horizontal 1 km long and 20 m wide lead was positioned in the center of the 5.12 km×5.12 km domain. In a second experiment, again with Δx = 20 m, the same lead was positioned in a domain twice this size; that is the boundaries are much farther from the lead in the second experiment. For both experiments, h_l is again equal to 2 m. The pressure fields around the lead are very similar (not shown) in the two experiments, with maximum pressures in the domain equal, respectively, to 36.22 and 36.36 kN m⁻¹ (a difference of ∼ 0.4 %). To avoid these boundary effects, we tend to position the important features in the center of the domain for the numerical experiments. For a numerical experiment to be valid, we require the simulated pressure in the first grid cells around the domain to be within 10 % of the pressure value specified at the boundaries.

6 Results

To limit the number of parameters that can be varied in the numerical experiments, the thickness of the level ice h_l is always set to 2 m. Furthermore, for all the experiments except the ones for the last figure, the normal stresses at the boundaries are always equal to −10 kN m⁻¹, while the shear stresses are set to 0. In other words, ${\mathit{\sigma }}^{\mathrm{W}}={\mathit{\sigma }}^{\mathrm{E}}={\mathit{\sigma }}^{\mathrm{S}}={\mathit{\sigma }}^{\mathrm{N}}$ = −10 kN m⁻¹ and ${\mathit{\tau }}^{\mathrm{W}}={\mathit{\tau }}^{\mathrm{E}}={\mathit{\tau }}^{\mathrm{S}}={\mathit{\tau }}^{\mathrm{N}}=\mathrm{0}\mathrm{kN}{\mathrm{m}}^{-\mathrm{1}}$.

6.1 Idealized sea ice experiments

Figure 4 Pressure (kN m⁻¹) field for a thickness field of 2 m everywhere except a 1 km long, 10 m wide horizontal lead in the middle of the domain (a; referred to as “single lead”). Pressure (kN m⁻¹) field for a thickness field of 2 m everywhere except a 1 km long, 10 m wide horizontal lead in the middle of the domain, a diagonal refrozen lead (h = 0.5 m), a smaller lead in the northwestern part of the domain and a 1 km ridge (h = 5 m in center, 2.5 m on each side) in the southeastern part of the domain (b, referred to as “many features”). For both experiments Δx = 10 m, and the normal stresses at the boundaries are −10 kN m⁻¹. c PDFs of the pressure for the “single lead” experiment (cyan) and the “many features” experiment (dashed magenta).

In a first set of experiments, we conduct idealized experiments to investigate the impact of sea ice features (leads, ridges, etc.) on the small-scale pressure field and especially on the maximum pressure. These experiments will give us insights and guide us for the second series of experiments with the idealized ship (see Sect. 6.2). Figure 4a shows the pressure field for a uniform sea ice cover with h_l=2 m except the presence of a 1 km long and 10 m wide lead. Large pressure is observed at the tips of the lead. In a second experiment, the same sea ice thickness conditions are used except that a smaller lead, a refrozen lead (with h=0.5 m) and a thick sea ice ridge (with h=5 m) are also positioned in the 5 km × 5 km domain. The pressure field for this latter experiment is shown in Fig. 4b. Figure 4c compares the PDFs of pressure for these two experiments. Looking at the PDFs and comparing Fig. 4a and b, one can notice that the other features are not associated with such high pressure values and that the maximum pressure is associated with the long 1 km lead. To further support this conclusion, note that the maximum pressure in the 5 km × 5 km domain is 42.57 kN m⁻¹ in the first experiment, while it is 42.59 kN m⁻¹ in the second one. In other words, the other features do not change our analysis; what really matters is the longest lead as it is in the vicinity of the longest lead that the largest stress concentration is found.

Figure 5 Maximum value of the shear stress invariant (a; kN m⁻¹) and of the pressure (b; kN m⁻¹) in the domain as a function of lead length for different values of the parameter P^* (P^* = 20 kN m⁻²: blue; P^* = 27.5 kN m⁻²: orange; P^* = 35 kN m⁻²: magenta). The thickness field is 2 m everywhere except the 10 m wide horizontal lead in the middle of the domain. The normal stresses at the boundaries are −10 kN m⁻¹.

Our results above suggest that only the longest lead needs to be considered for estimating the largest small-scale pressure. For a given h_l and stresses applied at the boundaries, there is more and more stress concentration when increasing the length of a lead. This is shown in Fig. 5 for three values of the parameter P^*. For short leads, the ice around the lead is able to sustain the stresses (the ice is rigid, that is in the viscous regime). This is why the three curves are very similar in Fig. 5a and b for short leads. However, for longer leads, there is more and more stress concentration. Some points of the ice close to the tips of the lead fail (i.e. the state of stress reaches the yield curve).

As the whole yield curve scales with the value of P^*, a larger P^* leads to larger maximum pressure and shear values. When increasing P^*, the maximum shear stress approaches asymptotically the shear strength (dashed lines in Fig. 5a; $e^{-\mathrm{1}}{h}_{\mathrm{l}}{P}^{*}/\mathrm{2}$). This asymptotic behaviour is less obvious for the pressure (Fig. 5b) as it is still far from the compressive strength (h_lP^*). A similar behaviour is observed when varying the ellipse aspect ratio (which modifies the shear strength). A smaller value of e leads to larger pressure values and larger shear stresses values (with a similar asymptotic behaviour) for long leads (not shown).

Figure 6 Pressure field with P^* = 27.5 kN m⁻² minus the pressure field with P^* = 20.0 kN m⁻² (kN m⁻¹). For both experiments, the thickness field is 2 m everywhere except a 1 km long, 10 m wide horizontal lead in the middle of the domain. The normal stresses at the boundaries are −10 kN m⁻¹.

While the average pressure in the domain is the same (10 kN m⁻¹) for all the values of P^*, the maximum pressure is enhanced as P^* increases (as shown in Fig. 5b). Comparing the pressure fields with P^* = 27.5 kN m⁻² and P^* = 20 kN m⁻² (see Fig. 6) for the same lead shows that the pressure fields around the lead are different over hundreds of metres. Moreover, the largest difference in the pressure fields are found at the tips of the lead; the pressure is much larger with P^* = 27.5 kN m⁻² than with P^* = 20 kN m⁻² in the vicinity of the tips.

Figure 7 Schematic of the coarse-graining procedure. The thickness field is defined at 10 m resolution (blue cells on the left). The thickness field at 20 m resolution is obtained by averaging the h values of the four 10 m cells contained in the 20 m one (purple cell). This procedure is repeated for the other lower resolutions. The same method is applied for the concentration A. The indices l,m are for the 10 m grid, while the indices $l^{*},m^{*}$ are for the 20 m grid.

We also investigate the evolution of the small-scale pressure field as a function of resolution. The h and A fields are defined at 10 m resolution. These fields h and A are, respectively, set to 2 m and 1.0 everywhere except for a 1 km long, 10 m wide lead in the middle of the domain with h=0 m and A=0. The model is run at resolutions of 10, 20, 40, … 1280 m. For these lower resolutions, the h and A fields are obtained through a coarse-graining procedure (see Fig. 7 for details).

Figure 8 All the values of pressure (a) and of the shear stress invariant (b) in the 5 km × 5 km domain as a function of resolution. The thickness field is 2 m everywhere except a 1 km long, 10 m wide horizontal lead in the middle of the domain. The thickness and concentration fields at the other resolutions are obtained through a coarse-graining procedure. The normal stresses at the boundaries are −10 kN m⁻¹.

All the values of p and q in the 5 km × 5 km domain are plotted as a function of Δx in Fig. 8a and b. The distribution of these small-scale stresses are non-symmetric (they are limited by 0 on one side) and are skewed toward large values. These results constitute another validation of the numerical framework as the distribution reduces to a single point equal to the large-scale values prescribed at the boundaries as Δx tends toward the horizontal dimension of the domain.

6.2 Experiments with an idealized ship

Figure 9 Pressure field at 10 m resolution when including a digitized ship 90 m long and 30 m wide (in grey). The thickness field is 2 m everywhere except a 600 m long lead behind the ship. The normal stresses at the boundaries are −10 kN m⁻¹. The whole domain is shown in (a), while (b) shows a zoom of the pressure field around the ship.

In a second set of experiments, we investigate the small-scale pressure field in the vicinity of a ship in heavy sea ice conditions and under compressive stresses. Importantly, we estimate the maximum pressure applied on the ship in different idealized experiments. The small-scale pressure field around a ship 90 m long and 30 m wide is investigated. We assume that the ship was navigating in level ice 2 m thick and that it is now beset. First, it is assumed that a lead (i.e. a channel) created by the ship is still open behind it over a distance of 600 m, while farther away the lead has been closed due to resulting sea ice convergence. The pressure field for this experiment is shown in Fig. 9a and with more details in Fig. 9b. Similar to our previous results without a ship, larger pressures are found at the tips of the lead. In fact, there is very large pressure on both sides of the ship, especially at the back of the ship. Numerical simulations of ships navigating in sea ice show larger pressure at the front of the ship (e.g. Kubat et al., 2010; Sayed et al., 2017). However, our results show the opposite for a ship that is beset. These results also suggest that by navigating in these compact ice conditions, the ship has generated these high-pressure conditions by creating a lead in its wake.

Figure 10 Maximum pressure (kN m⁻¹) on the ship as a function of the length of the lead behind the ship. The thickness field is 2 m everywhere except in the lead behind the ship; the thickness of the refrozen lead (h_rf) is 0 cm for the blue curve, it is 10 cm for the orange one, and it is 20 cm for the magenta one. The digitized ship is 90 m long and 30 m wide. The normal stresses at the boundaries are −10 kN m⁻¹.

A crucial aspect to consider here is the length of the lead behind the ship. Assuming the lead closes at a shorter distance from the ship should imply smaller pressure values (for the same pressure applied at the boundaries). This is indeed the case as demonstrated by the sensitivity study shown in Fig. 10 (blue curve). We also consider the case of a lead partially consolidated. In fact, we assume that the concentration (A) is 0 just behind the ship and that it increases linearly to 1.0 for a certain lead length. The (mean) thickness h of the ice is set equal to Ah_l. This appears to have a very small effect on our results (not shown) compared to the case with A=0 everywhere in the lead (blue curve in Fig. 10). This is due to the fact that the ice strength (see Eq. 7) decreases rapidly as A diminishes. However, if we consider that the ice in the lead is consolidating through thermodynamical growth (i.e. we set h to a small value in the lead behind the ship), we find that this notably reduces the stress applied on the ship. This can be seen with the orange and magenta curves in Fig. 10, which, respectively, correspond to thicknesses of 0.1 and 0.2 m for the refrozen lead.

Figure 11 Maximum pressure (kN m⁻¹) on the ship as a function of the pressure (p_b) prescribed at the boundaries. For both curves, it is assumed that the length of the lead (L) is 600 m for p_b = 0 kN m⁻¹ and that it decreases linearly as p_b increases. For the blue, the lead behind the ship is closed when p_b reaches 20 kN m⁻¹, while the lead is closed when p_b = 15 kN m⁻¹ for the magenta curve. The thickness field is 2 m everywhere except in the lead behind the ship, where the thickness is 0 m. The digitized ship is 90 m long and 30 m wide.

Figure 10 shows that, for a certain large-scale pressure applied at the boundaries, the length of the lead behind the ship has a strong impact on the maximum pressure applied on the ship. Even though it is unclear how long the lead should be for a given h_l and for a given large-scale pressure, it is realistic to suppose that the higher the pressure is at the boundaries, the shorter the lead is (i.e. it has closed over a certain distance due to the compressive stresses). Note that this is what is usually assumed for ships navigating in sea ice under compressive stresses (see for example Suominen and Kujala, 2012). In this last experiment, with results shown in Fig. 11, it is therefore assumed that the lead is 600 m long when the pressure at the boundaries is 0 kN m⁻¹ and that it decreases linearly to 0 m when the pressure reaches 20 kN m⁻¹ (blue curve) or when it reaches 15 kN m⁻¹ (magenta curve). The relation between the lead length L and the pressure at the boundaries p_b is therefore L = $-\mathrm{30}{p}_{\mathrm{b}}+\mathrm{600}$ for the blue curve and L = $-\mathrm{40}{p}_{\mathrm{b}}+\mathrm{600}$ for the magenta one. We therefore consider here that the lead has consolidated mechanically and that there is no thermodynamical growth. Figure 11 roughly exhibits three different regimes. In the first regime, for small pressure at the boundaries (i.e. long lead length), the maximum pressure on the ship is linearly related to p_b because it is independent of lead length. In the second regime, for large pressure at the boundaries (i.e. small lead length), the maximum pressure is most sensitive to the lead length, and we see again a linear dependence, with a negative slope, on p_b. In between, in the third regime, two compensating effects are playing out: a larger pressure at the boundaries causes the lead to be shorter, which decreases the stress concentration in the vicinity of the ship, making the maximum pressure weakly sensitive to the pressure at the boundary.

7 Conclusions

We have investigated how sea ice pressure could be downscaled at scales relevant for navigation. The distribution of pressure at small scales is associated with non-uniform sea ice conditions. The PDF of the small-scale pressure is non-symmetric (it is limited by 0 on one side) and is skewed toward large values. Our results indicate that what really determines the largest values of pressure is associated with defects, that is long leads. Because a lead itself is not able to sustain any stress (unless it has refrozen), the load is taken by the ice around the lead with especially large values of the stresses in the vicinity of the tips. A sensitivity study indicates that the small-scale distribution and maximum pressure are notably affected by the choice of the shear strength (e) and compressive strength (P^*) for the elliptical yield curve. This suggests that a different yield curve and different mechanical strength properties would also lead to significantly different results.

Idealized experiments with a digitized ship beset in heavy sea ice conditions show that stress concentration also occurs in the vicinity of the ship. In fact, our simulations show that the largest pressure applied on the ship is found on both sides at the back of the ship. These results are different than the ones of Kubat et al. (2010) and Sayed et al. (2017) because our idealized ship is beset, while they considered a digitized ship progressing in looser ice conditions.

We also argue that the ship itself is responsible for the strong concentration of stress on its side; the lead (or channel) it created by navigating in sea ice causes these large values of the stresses. Moreover, it is found that even a short lead causes pressure values notably larger than the pressure applied at the domain boundaries. The stresses on the ship should decrease as the ice in the lead consolidates (by either thermodynamical growth or closing of the lead). These conclusions highlight the difficulty of providing sub-grid-scale pressure forecasts for navigation applications as the important parameters (i.e. the length of the lead and the thickness of the refrozen ice) are not well constrained.

A significant advantage of our numerical framework is that stresses can be specified at the boundaries. However, it is also important to note its limitations. First, it can only calculate the pressure field for a ship beset in heavy sea ice conditions; it cannot simulate the sea ice stresses applied on a ship navigating in ice-infested waters (as in Kubat et al., 2010). Also, in reality, sea ice convergence can cause ridging, which can locally increase the yield strength of the ice. This strain hardening process was not considered in our numerical experiments; the maximum possible pressure in the domain is equal to P^*h_l. Another possible limitation of our numerical framework is that the ice is modelled as a continuum material rather than a collection of discrete particles. It would be very interesting to still apply stresses at the boundaries but to model the interactions between the sea ice and the idealized ship with a model based on discrete floes (e.g. Daley et al., 2014; Metrikin and Løset, 2013).

In our numerical experiments, the digitized ship is simply represented as a rigid body with no outflow and no slip boundary conditions applied on the contour. A more realistic numerical framework should also involve a better representation of ship–ice interactions. For example, as done by Kubat et al. (2010), a Coulomb friction condition could be applied on the ship contour.

Although the convergence criterion for the steady-state solution of the velocity field has been reached in all the numerical experiments described in this paper, it is worth mentioning that this came with tremendous difficulties for the JFNK solver; the non-linear convergence was really slow, and the solver failed on some occasions to reach the required drop in the Euclidean norm of the residual within the allowed 500 non-linear iterations. Lemieux et al. (2010) have already shown that the JFNK solver exhibits robustness issues as the grid is refined. In fact, it is really the number of unknowns that impacts the non-linear convergence and robustness of the JFNK solver. This clearly indicates that innovations and more sophisticated numerical methods (e.g. Mehlmann and Richter, 2017) would be very beneficial for the sea ice modelling community.

Appendix A: Stability analysis

A few observations were made concerning the numerical stability of our new numerical framework with stresses applied at the boundaries. In this appendix, we discuss and provide explanations for these limitations.

1. We have noticed that for a simulation to be numerically stable, σ^W should be equal to σ^E, σ^S should be equal to σ^N, and all the shear stresses at the boundaries should have the same value (i.e. ${\mathit{\tau }}^{\mathrm{W}}={\mathit{\tau }}^{\mathrm{E}}={\mathit{\tau }}^{\mathrm{S}}={\mathit{\tau }}^{\mathrm{N}}$). This can be easily understood by considering the ice in the domain as a single piece of ice. Assuming there is no shear stress, the sum of the forces applied on the ice along the x axis are

   $$\begin{array}{}\text{(A1)}& \sum F_x={\mathit{\sigma }}^{\mathrm{W}}\mathrm{\Delta }x-{\mathit{\sigma }}^{\mathrm{E}}\mathrm{\Delta }x.\end{array}$$

   For stability, ∑F_x should be 0 so that the ice does not accelerate indefinitely. This means that σ^W should be equal to σ^E. The same conclusion applies for σ^S and σ^N. Finally, a similar argument can be made for the shear stresses in terms of conservation of angular momentum. Interestingly, these conditions are the same ones found for the Cauchy tensor for the stresses at a point.

2. Dukowicz (1997) mentions that for numerical stability, the internal stresses should be 0 at open boundaries, while our simulations show that it is possible to obtain stable solutions with non-zero stresses prescribed at the boundaries. To understand this, we revisit the stability analysis described in Dukowicz (1997). As Dukowicz (1997), we consider a simplified 1D momentum equation. However, we also take into account the replacement method. With these considerations, our 1D momentum equation is given by

   $$\begin{array}{}\text{(A2)}& \mathit{\rho }h{\displaystyle \frac{\partial u}{\partial t}}={\displaystyle \frac{\partial \mathit{\sigma }}{\partial x}}.\end{array}$$

   For stability, the rheology term should dissipate kinetic energy (KE). To investigate this, we multiply Eq. (A2) by u and integrate it over the whole domain (x=0, i.e. the west side, and x=L, i.e. the east side of our domain).

   $$\begin{array}{}\text{(A3)}& \underset{\mathrm{0}}{\overset{L}{\int }}u\mathit{\rho }h{\displaystyle \frac{\partial u}{\partial t}}\mathrm{d}x=\underset{\mathrm{0}}{\overset{L}{\int }}u{\displaystyle \frac{\partial \mathit{\sigma }}{\partial x}}\mathrm{d}x\end{array}$$

   As advection and thermodynamics are not considered, the thickness field is constant in time, and we can write

   $$\begin{array}{}\text{(A4)}& \underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\mathit{\rho }h{u}^{\mathrm{2}}}{\mathrm{2}}}\right)\mathrm{d}x=\underset{\mathrm{0}}{\overset{L}{\int }}u{\displaystyle \frac{\partial \mathit{\sigma }}{\partial x}}\mathrm{d}x.\end{array}$$

   In 1D, $\mathit{\sigma }={\mathit{\alpha }}^{\mathrm{2}}\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{I}}-\mathit{\zeta }\mathrm{\Delta }$ with $\mathit{\zeta }=\frac{P_{\mathrm{p}}}{\mathrm{2}{\mathrm{\Delta }}^{*}}$, ${\mathrm{\Delta }}^{*}=min(\mathrm{\Delta },{\mathrm{\Delta }}_{\text{min}})$, ${\dot{\mathit{ϵ}}}_{\mathrm{I}}=\frac{\partial u}{\partial x}$ and $\mathrm{\Delta }=\mathit{\alpha }\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|$ with $\mathit{\alpha }=\sqrt{\mathrm{1}+e^{-\mathrm{2}}}$.

   The term on the right can be integrated by parts, that is

   $$\begin{array}{}\text{(A5)}& {\displaystyle }\underset{\mathrm{0}}{\overset{L}{\int }}u{\displaystyle \frac{\partial \mathit{\sigma }}{\partial x}}\mathrm{d}x={\left[u\mathit{\sigma }\right]}_{\mathrm{0}}^L-\underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\partial u}{\partial x}}\mathit{\sigma }\mathrm{d}x,\text{(A6)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\underset{\mathrm{0}}{\overset{L}{\int }}\left({\displaystyle \frac{\mathit{\rho }h{u}^{\mathrm{2}}}{\mathrm{2}}}\right)\mathrm{d}x& =u_L{\mathit{\sigma }}_L-u_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}\\ & -\underset{\mathrm{0}}{\overset{L}{\int }}\left({\mathit{\alpha }}^{\mathrm{2}}\mathit{\zeta }{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}-{\dot{\mathit{ϵ}}}_{\mathrm{I}}\mathit{\zeta }\mathrm{\Delta }\right)\mathrm{d}x,\end{array}\end{array}$$

   where the time derivative has been moved outside the integral because the region of integration is fixed (Dukowicz, 1997). Note that $u_L=u{|}_{x=L}$ (same idea for the other terms). The integral on the left is the total KE. From our results above we know that σ_L has to be equal to σ₀. By symmetry, we can also assume that $u_L=-u_{\mathrm{0}}$. Hence, with the definition of the viscous coefficient, we can then write Eq. (A6) as

   $$\begin{array}{}\text{(A7)}& {\displaystyle \frac{\partial }{\partial t}}KE=-\mathrm{2}{u}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{0}}-\underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{\mathit{\alpha }{P}_{\mathrm{p}}}{\mathrm{2}{\mathrm{\Delta }}^{*}}}\left(\mathit{\alpha }{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}-{\dot{\mathit{ϵ}}}_{\mathrm{I}}\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|\right)\mathrm{d}x.\end{array}$$

   For the second term on the right, $\left(\mathit{\alpha }{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}-{\dot{\mathit{ϵ}}}_{\mathrm{I}}\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|\right)={\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}(\mathit{\alpha }-\mathrm{1})$ if ${\dot{\mathit{ϵ}}}_{\mathrm{I}}$ is positive (divergence), while it is equal to ${\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}(\mathit{\alpha }+\mathrm{1})$ if ${\dot{\mathit{ϵ}}}_{\mathrm{I}}$ is negative (convergence). As α≥1, this means that the integral is always positive, and the term therefore always dissipates KE because of the minus sign in front of it. As opposed to the derivation of Dukowicz (1997), the replacement method is also considered here. Nevertheless, consistent with his results, we find that the second term on the right always dissipates KE.

   The stability therefore depends on the boundary term −2u₀σ₀. The worst condition happens when there are high compressive stresses at the boundaries. In this case, ${\mathit{\sigma }}_{\mathrm{0}}=-\left|{\mathit{\sigma }}_{\mathrm{0}}\right|<\mathrm{0}$ and u₀>0 such that $\mathrm{2}{u}_{\mathrm{0}}\left|{\mathit{\sigma }}_{\mathrm{0}}\right|$ is a source of KE. For high compressive stresses, we assume that the ice is in the plastic regime. To be able to evaluate the integral on the right in Eq. (A7), we also look at a simple case with P_p that is constant over the whole domain. With these assumptions we find

   $$\begin{array}{}\text{(A8)}& {\displaystyle \frac{\partial }{\partial t}}KE=\mathrm{2}{u}_{\mathrm{0}}\left|{\mathit{\sigma }}_{\mathrm{0}}\right|-{\displaystyle \frac{\mathit{\alpha }{P}_{\mathrm{p}}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{\displaystyle \frac{{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}}{\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|}}\mathrm{d}x+{\displaystyle \frac{P_{\mathrm{p}}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{\dot{\mathit{ϵ}}}_{\mathrm{I}}\mathrm{d}x.\end{array}$$

   With ${\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}/\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|=\left|{\dot{\mathit{ϵ}}}_{\mathrm{I}}\right|=-{\dot{\mathit{ϵ}}}_{\mathrm{I}}$, because ${\dot{\mathit{ϵ}}}_{\mathrm{I}}<\mathrm{0}$, we can then write

   $$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\partial }{\partial t}}KE=\mathrm{2}{u}_{\mathrm{0}}\left|{\mathit{\sigma }}_{\mathrm{0}}\right|+{\displaystyle \frac{\mathit{\alpha }{P}_{\mathrm{p}}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{\dot{\mathit{ϵ}}}_{\mathrm{I}}\mathrm{d}x+{\displaystyle \frac{P_{\mathrm{p}}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{\dot{\mathit{ϵ}}}_{\mathrm{I}}\mathrm{d}x.\end{array}$$

   With ${\int }_{\mathrm{0}}^L{\dot{\mathit{ϵ}}}_{\mathrm{I}}\mathrm{d}x={\int }_{\mathrm{0}}^L\frac{\partial u}{\partial x}\mathrm{d}x=u_L-u_{\mathrm{0}}=-\mathrm{2}{u}_{\mathrm{0}}$ we obtain

   $$\begin{array}{}\text{(A10)}& {\displaystyle \frac{\partial }{\partial t}}KE=\mathrm{2}{u}_{\mathrm{0}}\left|{\mathit{\sigma }}_{\mathrm{0}}\right|-(\mathit{\alpha }+\mathrm{1})P_{\mathrm{p}}{u}_{\mathrm{0}}.\end{array}$$

   This means that $\left|{\mathit{\sigma }}_{\mathrm{0}}\right|$ should be smaller than the compressive strength $(\mathit{\alpha }+\mathrm{1})P_{\mathrm{p}}/\mathrm{2}$ for the solution to be stable (i.e. the rheology term dissipates KE). A similar analysis can be conducted if we assume a tensile stress at the boundaries. In this case, we find that the stress $\left|{\mathit{\sigma }}_{\mathrm{0}}\right|$ at the boundaries should be smaller than the tensile strength $(\mathit{\alpha }-\mathrm{1})P_{\mathrm{p}}/\mathrm{2}$.

   To ensure numerical stability, Dukowicz (1997) mentions that the stresses should be 0 at the open boundaries. This is a stricter condition than the one we find here. We have indeed demonstrated that the solution is stable as long as the stresses prescribed at the boundaries are between the compressive and tensile strengths of the ice. Numerical experiments (in 2D) confirm this finding. For example, when prescribing normal stresses of −10 kN m⁻¹ on a uniform sea ice cover, the solution is stable if $h_{\mathrm{l}}>\mathrm{10}\mathrm{kN}{\mathrm{m}}^{-\mathrm{1}}/P^{*}$ (not shown).