Viscous-plastic Sea Ice Model Research Articles

On the sensitivity of sea ice deformation statistics to plastic damage

Abstract. We implement a plastic damage parametrization, distinct from the elastic damage in the elasto-brittle framework, in the standard viscous–plastic (VP) sea ice model to disentangle its effect from resolved model physics (visco-plastic with and without damage) on its ability to reproduce observed scaling laws of deformation. To this end, we compare scaling properties and multifractality of simulated divergence and shear strain rate, as proposed in the Sea Ice Rheology Experiment (SIREx) studies, with those derived from the RADARSAT Geophysical Processor System (RGPS). Results show that including a plastic damage parametrization in the standard viscous–plastic model increases the spatial but decreases the temporal localization of simulated linear kinematic features (LKFs) and brings all spatial deformation rate statistics in line with observations from RGPS without the need to increase the mechanical shear strength of sea ice as recently proposed for lower-resolution viscous–plastic sea ice models. In fact, including damage with a healing timescale of th=30 d and an increased mechanical strength unveils multifractal behavior that does not fit the theory. Therefore, a plastic damage parametrization is a powerful tuning knob affecting the deformation statistics of viscous–plastic sea ice.

A dynamical core based on a discontinuous Galerkin method for higher-order finite-element sea ice modeling

Abstract. The ability of numerical sea ice models to reproduce localized deformation features associated with fracture processes is key for an accurate representation of the ice dynamics and of dynamically coupled physical processes in the Arctic and Antarctic. Equally key is the capacity of these models to minimize the numerical diffusion stemming from the advection of these features to ensure that the associated strong gradients persist in time, without the need to unphysically re-inject energy for re-localization. To control diffusion and improve the approximation quality, we present a new numerical core for the dynamics of sea ice that is based on higher-order finite-element discretizations for the momentum equation and higher-order discontinuous Galerkin methods for the advection. The mathematical properties of this core are discussed, and a detailed description of an efficient shared-memory parallel implementation is given. In addition, we present different numerical tests and apply the new framework to a benchmark problem to quantify the advantages of the higher-order discretization. These tests are based on Hibler's viscous–plastic sea ice model, but the implementation of the developed framework in the context of other physical models reproducing a strong localization of the deformation is possible.

Toward a method for downscaling sea ice pressure for navigation purposes

Abstract. 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.

On the calculation of normalized viscous–plastic sea ice stresses

Abstract. Calculating and plotting the normalized states of stress for viscous–plastic sea ice models is a common diagnostic for evaluating the numerical convergence and the physical consistency of a numerical solution. Researchers, however, usually do not explain how they calculate the normalized stresses. Here, we argue that care must be taken when calculating and plotting the normalized states of stress. A physically consistent and numerically converged solution should exhibit normalized stresses that are inside (viscous) or on (plastic) the normalized yield curve. To do so, two possible mistakes need to be avoided. First, when using an implicit solver, normalized stresses should be computed from viscous coefficients and replacement pressure calculated using the previous numerical iterate and the strain rates at the numerator calculated from the latest iterate. Calculating the stresses only from the latest iterate falsely indicates that the solution has numerically converged. Second, for both implicit and explicit (i.e., the EVP) solvers, the stresses should be normalized by the ice strength and not by the replacement pressure. Using the latter, normalized states of stress only lie on the yield curve (i.e., falsely indicating there are no viscous states of stress).

The dependence of energy dissipation on spatial resolution in a viscous-plastic sea-ice model

We present sea-ice kinetic energy budgets to quantify the relative importance of the various energy sinks in a viscous-plastic sea-ice model. To this end, we study two idealized model domains where energy dissipation associated with shear and axial (ridge/lead building) deformation can be analyzed independently. We find that when only shear deformation is present - either at the domain boundary induced by the no-slip boundary condition or within the model domain induced by gradients in the surface air stress - the energy dissipated through friction reduces in relative importance as the spatial resolution of the model is increased. In the limit where the spatial resolution tends to zero, the simulated sea ice drift tends to the analytical solution - giving us confidence in the numerical implementation of the governing differential equation. Increasing spatial resolution leads to a localization of deformation along the shear lines effectively increasing the area over which energy is input by the wind which is not compensated for by frictional shear dissipation. For instance at 40 km spatial resolution, 64% and 29% of the input power is dissipated through shear deformation and water drag respectively, while at 5 km spatial resolution 54% and 43% of the input power is dissipated by the respective processes. These values approach the respective values of 53% and 47% found analytically for this particular model configuration. The overall result is a 64% increase in the domain total sea-ice kinetic energy when the spatial resolution is increased from 40 km to 5 km due to the finer representation of shear lines. In convergence, the mean kinetic energy and potential energy do not depend meaningfully on the spatial resolution. In this case, the structure of the thickness and concentration fields effectively sets the velocity gradient near the boundary provided that the plastic deformation wave associated with the ridge building process is resolved.

A modified global Newton solver for viscous-plastic sea ice models

We present and analyze a modified Newton solver, the so called operator-related damped Jacobian method, with a line search globalization for the solution of the strongly nonlinear momentum equation in a viscous-plastic (VP) sea ice model.Due to large variations in the viscosities, the resulting nonlinear problem is very difficult to solve. The development of fast, robust and converging solvers is subject to present research. There are mainly three approaches for solving the nonlinear momentum equation of the VP model, a fixed-point method denoted as Picard solver, an inexact Newton method and a subcycling procedure based on an elastic-viscous-plastic model approximation. All methods tend to have problems on fine meshes by sharp structures in the solution. Convergence rates deteriorate such that either too many iterations are required to reach sufficient accuracy or convergence is not obtained at all.To improve robustness globalization and acceleration approaches, which increase the area of fast convergence, are needed. We develop an implicit scheme with improved convergence properties by combining an inexact Newton method with a Picard solver. We derive the full Jacobian of the viscous-plastic sea ice momentum equation and show that the Jacobian is a positive definite matrix, guaranteeing global convergence of a properly damped Newton iteration. We compare our modified Newton solver with line search damping to an inexact Newton method with established globalization and acceleration techniques. We present a test case that shows improved robustness of our new approach, in particular on fine meshes.

On Well-Posedness of the Viscous–Plastic Sea Ice Model

Abstract This study examines the well-posedness of the initial-value problems that arise in common models of sea ice. The model equations describe the balance of linear momentum combined with simplified thermodynamics represented by two continuity equations for effective ice thickness and ice concentration. The constitutive model for sea ice is given by two possible variants of the viscous–plastic model: the viscous–plastic model with pressure replacement and the viscous–plastic model with pressure replacement plus a tensile cutoff. The authors identify regimes of well- and ill-posedness for both models in one and two space dimensions. In one space dimension, the study finds that the viscous–plastic model and viscous–plastic model with pressure replacement behave similarly: there is ill-posedness when the divergent flow rate is larger than a minimum value. On the other hand, the viscous–plastic model with pressure replacement plus a tensile cutoff is ill-posed for all divergent flows. In two space dimensions the analysis is inconclusive for the viscous–plastic model with pressure replacement, but with the tensile cutoff the problem is ill-posed for certain divergent flows. The authors also discuss energy bounds and the difference between ill-posedness and stability of a solution. The study shows by examples that boundedness of solutions does not imply well-posedness and that it is possible for well-posed problems to have unstable solutions. The analysis shows that previous arguments in the literature, which state that a bound on the energy in sea ice models provides control over ill-posedness, are flawed.

Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton–Krylov method

We have implemented the Jacobian-free Newton–Krylov (JFNK) method to solve the sea ice momentum equation with a viscous-plastic (VP) formulation. The JFNK method has many advantages: the system matrix (the Jacobian) does not need to be formed and stored, the method is parallelizable and the convergence can be nearly quadratic in the vicinity of the solution. The convergence rate of our JFNK implementation is characterized by two phases: an initial phase with slow convergence and a fast phase for which the residual norm decreases significantly from one Newton iteration to the next. Because of this fast phase, the computational gain of the JFNK method over the standard solver used in existing VP models increases with the required drop in the residual norm (termination criterion). The JFNK method is between 3 and 6.6 times faster (depending on the spatial resolution and termination criterion) than the standard solver using a preconditioned generalized minimum residual method. Resolutions tested in this study are 80, 40, 20 and 10 km. For a large required drop in the residual norm, both JFNK and standard solvers sometimes do not converge. The failure rate for both solvers increases as the grid is refined but stays relatively small (less than 2.3% of failures). With increasing spatial resolution, the velocity gradients (sea ice deformations) get more and more important. Nonlinear solvers such as the JFNK method tend to have difficulties when there are such sharp structures in the solution. This lack of robustness of both solvers is however a debatable problem as it mostly occurs for large required drops in the residual norm. Furthermore, when it occurs, it usually affects only a few grid cells, i.e., the residual is small for all the velocity components except in very localized regions. Globalization approaches for the JFNK solver, such as the line search method, have not yet proven to be successful. Further investigation is needed.

An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids

The elastic–viscous–plastic (EVP) sea ice rheology has been introduced in the large-scale Louvain-la-Neuve sea-Ice Model, version 2 (LIM2), and its performance has been evaluated. Centred difference versions of the rheology have been implemented on both an Arakawa B grid and a C grid, and their performance have been intercompared in coupled simulations with the Nucleus for European Modelling of the Ocean (NEMO) model. Integrations with both implementations lead to fairly similar results which compare well with observations and with previous LIM simulations. The C grid version, however, offers a number of advantages: (a) easier ice coupling with NEMO, which is itself defined on a C grid; (b) possibility of representing ice transport across one-cell-wide straits and passages; (c) better representation of inertial-plastic compressive waves. For these reasons, we recommend the use of the C grid EVP formulation of the ice dynamics in future LIM applications.

Modeling Linear Kinematic Features in Sea Ice

Abstract Sea ice deformation is localized in narrow zones of high strain rate that extend hundreds of kilometers, for example, across the Arctic Basin. This paper demonstrates that these failure zones may be modeled with a viscous–plastic sea ice model, using an isotropic rheology. If the ice is assumed to be heterogeneous at the grid scale, and allowed to weaken in time, intersecting failure zones propagate across the region. The direction of failure propagation depends upon the stress applied to the ice (wind stress and boundary conditions) and the rheological model describing plastic failure of the ice. The spacing between failure zones is controlled by the magnitude of the wind stress and the distribution describing spatial variability of ice strength. Sea ice motion and deformation oscillate at close to a 12-h period throughout the Arctic and Antarctic pack ice. This oscillation is found at all spatial scales from hundreds of kilometers to the lead scale. It is shown that with an inertial embedded model, sea ice deformation rotates between pairs of fault patterns with a semidiurnal period. It is well known that linear zones of deformation exist at many spatial scales throughout the Arctic Basin. The model presented in this paper may be scaled to simulate these features.

A strength implicit correction scheme for the viscous-plastic sea ice model

With increasing awareness of the role of the polar regions in global climate, efficient and accurate simulation of sea ice is an important issue. The numerical solution procedure for the viscous-plastic sea ice model is examined. Recent developments are drawn upon and new discretisation practices introduced that improve efficiency and the accuracy of the simulation. The main components of the new solution procedure are bounded and conservative finite volume discretisation and a reinterpretation of the model structure, including the introduction of a method to resolve momentum-ice strength coupling. It is found that traditional segregated momentum solution procedures do not ensure mass conservation, which may introduce substantial thickness errors in Arctic sea ice simulations. The strength implicit correction scheme ensures mass conservation for a similar computational expense as previous, non-conservative methods.

Parallel climate model (PCM) control and transient simulations

The Department of Energy (DOE) supported Parallel Climate Model (PCM) makes use of the NCAR Community Climate Model (CCM3) and Land Surface Model (LSM) for the atmospheric and land surface components, respectively, the DOE Los Alamos National Laboratory Parallel Ocean Program (POP) for the ocean component, and the Naval Postgraduate School sea-ice model. The PCM executes on several distributed and shared memory computer systems. The coupling method is similar to that used in the NCAR Climate System Model (CSM) in that a flux coupler ties the components together, with interpolations between the different grids of the component models. Flux adjustments are not used in the PCM. The ocean component has 2/3° average horizontal grid spacing with 32 vertical levels and a free surface that allows calculation of sea level changes. Near the equator, the grid spacing is approximately 1/2° in latitude to better capture the ocean equatorial dynamics. The North Pole is rotated over northern North America thus producing resolution smaller than 2/3° in the North Atlantic where the sinking part of the world conveyor circulation largely takes place. Because this ocean model component does not have a computational point at the North Pole, the Arctic Ocean circulation systems are more realistic and similar to the observed. The elastic viscous plastic sea ice model has a grid spacing of 27 km to represent small-scale features such as ice transport through the Canadian Archipelago and the East Greenland current region. Results from a 300 year present-day coupled climate control simulation are presented, as well as for a transient 1% per year compound CO2 increase experiment which shows a global warming of 1.27 °C for a 10 year average at the doubling point of CO2 and 2.89 °C at the quadrupling point. There is a gradual warming beyond the doubling and quadrupling points with CO2 held constant. Globally averaged sea level rise at the time of CO2 doubling is approximately 7 cm and at the time of quadrupling it is 23 cm. Some of the regional sea level changes are larger and reflect the adjustments in the temperature, salinity, internal ocean dynamics, surface heat flux, and wind stress on the ocean. A 0.5% per year CO2 increase experiment also was performed showing a global warming of 1.5 °C around the time of CO2 doubling and a similar warming pattern to the 1% CO2 per year increase experiment. El Nino and La Nina events in the tropical Pacific show approximately the observed frequency distribution and amplitude, which leads to near observed levels of variability on interannual time scales.

The role of plastic ice interaction in marginal ice zone dynamics

Under appropriate conditions, the nonlinear nature of plastic ice interaction together with a nonlinear coupling between ice thickness characteristics and ice rheology can substantially modify the character of marginal ice zone dynamics. This paper examines the steady state ramifications of these nonlinearities by using a one‐dimensional simplification of Hibler's (1979) two‐level viscous plastic sea ice model. A series of idealized small‐scale simulations (4‐km resolution) are carried out with the model formulated in a moving Lagrangian grid in order to remove diffusion effects. Analytic solutions for the equilibrium plastic adjustment case are also constructed. The results show that if the ice thickness distribution is allowed to equilibrate in response to a constant wind field, the thickness strength coupling will yield a sharp ice edge, with the compactness dropping rapidly to zero near the ice margin. Despite this compactness falloff, the steady adjusted ice velocity is nearly constant across the marginal ice zone. However, the adjustment to a steady state constant velocity can be preceded by a period of shearing, which will abate as the ice strengthens by ridging. With an increase in forcing near the edge, the ice velocity scales with the wind velocity but with a speed less than free drift. In the absence of water drag and wind stress the model exhibits a very slow outward creep but yields no ice edge jet effects. The ramifications of plastic interaction on upwelling and downwelling are also examined. With the same drag coefficients over ice and ocean, and for a north‐south ice edge, the model yields (in the northern hemisphere) upwelling for northward winds, downwelling for southward winds, and no effect for off‐ice winds.