On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations

This paper describes the sea ice component of the Massachusetts Institute of Technology general circulation model (MITgcm); it presents example Arctic and Antarctic results from a realistic, eddy-admitting, global ocean and sea ice configuration; and it compares B-grid and C-grid dynamic solvers and other numerical details of the parameterized dynamics and thermodynamics in a regional Arctic configuration. Ice mechanics follow a viscous-plastic rheology and the ice momentum equations are solved numerically using either line-successive-over-relaxation (LSOR) or elastic–viscous-plastic (EVP) dynamic models. Ice thermodynamics are represented using either a zero-heat-capacity formulation or a two-layer formulation that conserves enthalpy. The model includes prognostic variables for snow thickness and for sea ice salinity. The above sea ice model components were borrowed from current generation climate models but they were reformulated on an Arakawa C grid in order to match the MITgcm oceanic grid and they were modified in many ways to permit efficient and accurate automatic differentiation. Both stress tensor divergence and advective terms are discretized with the finite-volume method. The choice of the dynamic solver has a considerable effect on the solution; this effect can be larger than, for example, the choice of lateral boundary conditions, of ice rheology, and of ice-ocean stress coupling. The solutions obtained with different dynamic solvers typically differ by a few cm s−1 in ice drift speeds, 50 cm in ice thickness, and order 200 km3 yr−1 in freshwater (ice and snow) export out of the Arctic.