Jean-François Lemieux

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 10km. 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.

Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the sea‐ice momentum equation

[1] We introduce the preconditioned generalized minimum residual (GMRES) method, along with an outer loop (OL) iteration to solve the sea-ice momentum equation. The preconditioned GMRES method is the linear solver. GMRES together with the OL is used to solve the nonlinear momentum equation. The GMRES method has low storage requirements, and it is computationally efficient and parallelizable. It was found that the preconditioned GMRES method is about 16 times faster than a stand-alone successive overrelaxation (SOR) solver and three times faster than a stand-alone line SOR (LSOR). Unlike stand-alone SOR and stand-alone LSOR, the cpu time needed by the preconditioned GMRES method for convergence weakly depends on the relaxation parameter when it is smaller than the optimal value. Results also show that with a 6-hour time step, the free drift velocity field is a better initial guess than the previous time step solution. For GMRES, the symmetry of the system matrix is not a prerequisite. The Coriolis term and the off-diagonal part of the water drag term can then be treated implicitly. The implicit treatment eliminates an instability characterized by a residual oscillation in the total kinetic energy of the ice pack that can be present when these off-diagonal terms are handled explicitly. Treating these terms explicitly prevents one from obtaining a high-accuracy solution of the sea-ice momentum equation unless a corrector step is applied. In fact, even after a large number of OL iterations, errors in the drift of the same magnitude as the drift itself can be present when these terms are treated explicitly.

Intercomparison of the Arctic sea ice cover in global ocean–sea ice reanalyses from the ORA-IP project

AbstractOcean–sea ice reanalyses are crucial for assessing the variability and recent trends in the Arctic sea ice cover. This is especially true for sea ice volume, as long-term and large scale sea ice thickness observations are inexistent. Results from the Ocean ReAnalyses Intercomparison Project (ORA-IP) are presented, with a focus on Arctic sea ice fields reconstructed by state-of-the-art global ocean reanalyses. Differences between the various reanalyses are explored in terms of the effects of data assimilation, model physics and atmospheric forcing on properties of the sea ice cover, including concentration, thickness, velocity and snow. Amongst the 14 reanalyses studied here, 9 assimilate sea ice concentration, and none assimilate sea ice thickness data. The comparison reveals an overall agreement in the reconstructed concentration fields, mainly because of the constraints in surface temperature imposed by direct assimilation of ocean observations, prescribed or assimilated atmospheric forcing and assimilation of sea ice concentration. However, some spread still exists amongst the reanalyses, due to a variety of factors. In particular, a large spread in sea ice thickness is found within the ensemble of reanalyses, partially caused by the biases inherited from their sea ice model components. Biases are also affected by the assimilation of sea ice concentration and the treatment of sea ice thickness in the data assimilation process. An important outcome of this study is that the spatial distribution of ice volume varies widely between products, with no reanalysis standing out as clearly superior as compared to altimetry estimates. The ice thickness from systems without assimilation of sea ice concentration is not worse than that from systems constrained with sea ice observations. An evaluation of the sea ice velocity fields reveals that ice drifts too fast in most systems. As an ensemble, the ORA-IP reanalyses capture trends in Arctic sea ice area and extent relatively well. However, the ensemble can not be used to get a robust estimate of recent trends in the Arctic sea ice volume. Biases in the reanalyses certainly impact the simulated air–sea fluxes in the polar regions, and questions the suitability of current sea ice reanalyses to initialize seasonal forecasts.

A basal stress parameterization for modeling landfast ice

Current large-scale sea ice models represent very crudely or are unable to simulate the formation, maintenance and decay of coastal landfast ice. We present a simple landfast ice parameterization representing the effect of grounded ice keels. This parameterization is based on bathymetry data and the mean ice thickness in a grid cell. It is easy to implement and can be used for two-thickness and multithickness category models. Two free parameters are used to determine the critical thickness required for large ice keels to reach the bottom and to calculate the basal stress associated with the weight of the ridge above hydrostatic balance. A sensitivity study was conducted and demonstrates that the parameter associated with the critical thickness has the largest influence on the simulated landfast ice area. A 6 year (2001–2007) simulation with a 20 km resolution sea ice model was performed. The simulated landfast ice areas for regions off the coast of Siberia and for the Beaufort Sea were calculated and compared with data from the National Ice Center. With optimal parameters, the basal stress parameterization leads to a slightly shorter landfast ice season but overall provides a realistic seasonal cycle of the landfast ice area in the East Siberian, Laptev and Beaufort Seas. However, in the Kara Sea, where ice arches between islands are key to the stability of the landfast ice, the parameterization consistently leads to an underestimation of the landfast area.

The Granular Sea Ice Model in Spherical Coordinates and Its Application to a Global Climate Model

Abstract The granular sea ice model (GRAN) from Tremblay and Mysak is converted from Cartesian to spherical coordinates. In this conversion, the metric terms in the divergence of the deviatoric stress and in the strain rates are included. As an application, the GRAN is coupled to the global Earth System Climate Model from the University of Victoria. The sea ice model is validated against standard datasets. The sea ice volume and area exported through Fram Strait agree well with values obtained from in situ and satellite-derived estimates. The sea ice velocity in the interior Arctic agrees well with buoy drift data. The thermodynamic behavior of the sea ice model over a seasonal cycle at one location in the Beaufort Sea is validated against the Surface Heat Budget of the Arctic Ocean (SHEBA) datasets. The thermodynamic growth rate in the model is almost twice as large as the observed growth rate, and the melt rate is 25% lower than observed. The larger growth rate is due to thinner ice at the beginning of ...

Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges

In some coastal regions of the Arctic Ocean, grounded ice ridges contribute to stabilizing and maintaining a landfast ice cover. Recently, a grounding scheme representing this effect on sea ice dynamics was introduced and tested in a viscous-plastic sea ice model. This grounding scheme, based on a basal stress parameterization, improves the simulation of landfast ice in many regions such as in the East Siberian Sea, the Laptev Sea and along the coast of Alaska. Nevertheless, in some regions such as the Kara Sea, the area of landfast ice is systematically underestimated. This indicates that another mechanism such as ice arching is at play for maintaining the ice cover fast. To address this problem, the combination of the basal stress parameterization and tensile strength is investigated using a 0.25° pan-Arctic CICE-NEMO configuration. Both uniaxial and isotropic tensile strengths notably improve the simulation of landfast ice in the Kara Sea but also in the Laptev Sea. However, the simulated landfast ice season for the Kara Sea is too short compared to observations. This is especially obvious for the onset of the landfast ice season which systematically occurs later in the model and with a slower build up. This suggests that improvements to the sea ice thermodynamics could reduce these discrepancies with the data. This article is protected by copyright. All rights reserved.

The effects of plastic waves on the numerical convergence of the viscousplastic and elasticviscousplastic sea-ice models

The plastic wave speed is derived from the linearized 1-D version of the widely used viscousplastic (VP) and elasticviscousplastic (EVP) sea-ice models. CourantFriedrichsLewy (CFL) conditions are derived using the propagation speed of the wave. 1-D numerical experiments of the VP, EVP and EVP models successfully recreate a reference solution when the CFL conditions are satisfied, in agreement with the theory presented. The IMplicitEXplicit (IMEX) method is shown to effectively alleviate the plastic wave CFL constraint on the timestep in the implicitly solved VP model in both 1-D and 2-D. In 2-D, the EVP and EVP models show first order error in the simulated velocity field when the plastic wave is not resolved. EVP simulations are performed with various advective timestep, number of subcycles, and elastic-wave damping timescales. It is found that increasing the number of subcycles beyond that needed to resolve the elastic wave does not improve the quality of the solution. It is found that reducing the elastic wave damping timescale reduces the spatial extent of first order errors cause by the unresolved plastic wave. Reducing the advective timestep so that the plastic wave is resolved also reduces the velocity error in terms of magnitude and spatial extent. However, the parameter set required for convergence to within the error bars of satellite (RGPS) deformation fields is impractical for use in climate model simulations. The behavior of the EVP method is analogous to that of the EVP method except that it is not possible to reduce the damping timescale with =.

What historical landfast ice observations tell us about projected ice conditions in Arctic Archipelagoes and marginal seas under anthropogenic forcing

Arctic landfast ice extent and duration from observations, ice assimilations, ocean re-analyses and coupled models are examined. From observations and assimilations, it is shown that in areas where landfast ice conditions last more than 5 months the first-year ice grows typically to more than 2 m and is rarely less than 1 m. The observed spatial distribution of landfast ice closely matches assimilation products but less so for ocean re-analyses and coupled models. Although models generally struggle to represent the landfast ice necessary to emulate the observed import/export of sea ice in regions favourable to landfast ice conditions, some do exhibit both a realistic climatology and a realistic decline of landfast ice extent under an anthropogenic forcing scenario. In these more realistic simulations, projections show that an extensive landfast ice cover should remain for at least 5 months of the year well until the end of the 21st century. This is in stark contrast with the simulations that have an unrealistic emulation of landfast ice conditions. In these simulations, slow and packed ice condi tions shrink markedly over the same period. In all simulations and in areas with landast ice that last more than 5 months, the end-of-winter sea ice thickness remains between 1 m and 2 m well beyond the second half of the century. It is concluded that in the current generation of climate models, projections of winter sea ice conditions in the Canadian Arctic Archipelago and the Laptev Sea are overly sensitive to the representation of landfast ice conditions and that ongoing development in landfast ice parametrization will likely better constrain these projections.

Implementation of Newton's method with an analytical Jacobian to solve the 1D sea ice momentum equation

New numerical solvers are being considered in response to the rising computational cost of properly solving the sea ice momentum equation at high resolution. The Jacobian free version of Newtons method has allowed models to obtain the converged solution faster than other implicit solvers used previously. To further improve on this recent development, the analytical Jacobian of the 1D sea ice momentum equation is derived and used inside Newtons method. The results are promising in terms of computational efficiency. Although robustness remains an issue for some test cases, it is improved compared to the Jacobian free approach. In order to make use of the strong points of both the new and Jacobian free methods, a hybrid preconditioner using the Picard and Jacobian matrices to improve global and local convergence, respectively, is also introduced. This preconditioner combines the robustness and computational efficiency of the previously used preconditioning matrices when solving the sea ice momentum equation.