The Louvain-La-Neuve sea ice model LIM3.6: global and regional capabilities

Abstract

The new 3.6 version of the Louvain-la-Neuve sea ice model (LIM) is presented, as integrated in the most recent stable release of Nucleus for European Modelling of the Ocean (NEMO) (3.6). The release will be used for the next Climate Model Inter-comparison Project (CMIP6). Developments focussed around three axes: improvements of robustness, versatility and sophistication of the code, which involved numerous changes. Robustness was improved by enforcing exact conservation through the inspection of the different processes driving the air–ice–ocean exchanges of heat, mass and salt. Versatility was enhanced by implementing lateral boundary conditions for sea ice and more flexible ice thickness categories. The latter includes a more practical computation of category boundaries, parameterizations to use LIM3.6 with a single ice category and a flux redistributor for coupling with atmospheric models that cannot handle multiple sub-grid fluxes. Sophistication was upgraded by including the effect of ice and snow weight on the sea surface. We illustrate some of the new capabilities of the code in two standard simulations. One is an ORCA2-LIM3 global simulation at a nominal 2${}^{\circ }$ resolution, forced by reference atmospheric climatologies. The other one is a regional simulation at 2 km resolution around the Svalbard Archipelago in the Arctic Ocean, with open boundaries and tides. We show that the LIM3.6 forms a solid and flexible base for future scientific studies and model developments.

Introduction

Sea ice covers 3–6 % of the Earth's surface and is characterized by ample seasonal variations, making it one of the most influential geophysical features in the Earth system (Comiso, 2010). Mostly because of its high albedo and thermal insulation power, sea ice affects the weather and more generally the climate (e.g., Budkyko, 1969; Vihma, 2014). The seasonal cycle of ice growth and melt strongly impacts the vertical upper ocean density structure and is a key driver of the ocean circulation at a global scale through dense water formation (Aagaard and Carmack, 1989; Goosse and Fichefet, 1999). Sea ice also influences marine primary productivity and carbon export to depth (e.g. Thomas and Dieckmann, 2010; Vancoppenolle et al., 2013), and constitutes a habitat for Arctic and Antarctic fauna, including specific microbial, birds and mammal species (Croxall et al., 2002; Atkinson et al., 2004).

Given the difficulty to observe polar regions, numerical modelling is essential to understand sea ice processes and their influence on the other components of the Earth system. Indeed, a sea ice component is presently included in virtually all ocean and Earth modelling systems (e.g. Flato et al., 2013; Danabasoglu et al., 2014). The contemporary use of sea ice models encompasses a wide range of applications, from large-scale climate to small-scale process studies and operational forecasts. The physical processes at stake need to be well resolved at the appropriate spatial and temporal scales. Hence, sea ice models must be both physically reliable and versatile in a wide range of scales, at a reasonable computational cost (e.g. Hunke et al., 2010).

In order to match these constraints, a number of changes have been made into the Louvain-la-Neuve sea ice model (LIM3; Vancoppenolle et al., 2009a), leading to the 3.6 version of the code. Along with the interface for Community Ice CodE (CICE) (Hunke et al., 2013), LIM3.6 is now the reference sea ice model in the Nucleus for European Modelling of the Ocean (NEMO) in its 3.6 version just released in June 2015. NEMO-LIM3.6 is expected to have a long life time, as it will form the base of the ocean and sea ice representation in several forthcoming Earth system models for the Coupled Model Inter-comparison Project 6 (CMIP6) (Meehl et al., 2014): the Institut Pierre-et-Simon Laplace (IPSL) Earth system model (Dufresne et al., 2013), EC-Earth (Hazeleger et al., 2010) and Centro Euro-Mediterraneo sui Cambiamenti Climatici Climate Model (CMCC-CM) (Scoccimarro et al., 2011). Therefore, we found it timely and appropriate to present the new characteristics and possibilities made possible by LIM3.6 in this paper.

The modifications made to LIM mainly improve the robustness, versatility and sophistication of the code, aiming at satisfying the needs of a large community of users. A major goal was to reach an exact conservation of mass, heat and salt, which is essential for climate simulations but was not satisfied in LIM3 until now. For that purpose, the time stepping scheme was reshaped and several minor conservation leaks were found and corrected. New capabilities have also been developed: open-boundary conditions for sea ice (which enables regional studies in ice-covered areas), more flexible thickness category boundaries, mono-category parameterizations, more realistic ice–ocean interactions, and more flexible ice–atmosphere exchanges.

The representation of sea ice physics in LIM is described in Sect. 2. Section 3 is dedicated to the new developments of the sea ice model. Some of these developments are illustrated in two simulations using the latest stable release of NEMO-LIM: a large-scale global 2${}^{\circ }$-resolution configuration (Sect. 4); and a regional 2 km resolution configuration around the Svalbard Archipelago (Sect. 5), a region well-suited to study various sea ice regimes as well as transient events such as polynyas. Conclusions and perspectives are presented in Sect. 6.

Model description

LIM was originally a B-grid sea ice model developed by Fichefet and Morales-Maqueda (1997), including ice dynamics based on the viscous-plastic (VP) rheology (Hibler III, 1979), the three-layer thermodynamic formulation of Semtner Jr. (1976), the second-order moment-conserving advection scheme of Prather (1986) and various sea ice physical parameterizations. Some years later LIM became LIM2 when it was rewritten in Fortran 90 and coupled to Ocean Parallelise (OPA), a C-grid, finite difference, hydrostatic, primitive equation ocean general circulation model (Madec, 2008). LIM2 was later on integrated into the NEMO system, for the global reference configuration ORCA2-LIM (Timmermann et al., 2005).

Recently, LIM was improved towards a better account of sub-grid-scale physics, giving birth to LIM3 (Vancoppenolle et al., 2009a, b). LIM3, as other multi-category models (e.g. CICE; Hunke et al., 2013), is based on the Arctic Ice Dynamics Joint Experiment (AIDJEX) framework (Coon et al., 1974). LIM3 mostly differs from other multi-category models in terms of parameterizations and implementation details. The new features of LIM3 are mainly multiple ice categories to represent the sub-grid-scale ice thickness distribution (Thorndike et al., 1975), multi-layer halo-thermodynamics including brine dynamics and their impact on thermal properties and ice–ocean salt exchanges (Vancoppenolle et al., 2009b) and a C-grid formulation (Bouillon et al., 2009) for ice dynamics using the modified elastic-viscous-plastic (EVP) rheology (Bouillon et al., 2013), instead of the more computationally expensive VP rheology (Hibler III, 1979).

Conservation of area and ice thickness categories

To account for unresolved sub-grid-scale variations in ice thickness ($h)$, the state of sea ice is given by a thickness distribution function $g(x,y,h,t)$ (Thorndike et al., 1975), defined as the limit $$g=\underset{\text{d}h\to \mathrm{0}}{lim}\frac{\text{d}A}{\text{d}h},$$ where $\text{d}A$ is the areal fraction of a small control surface with thickness between $h$ and $h+\text{d}h$.

Invoking continuity, the conservation of area can be written as $$\frac{\partial g}{\partial t}=-\mathrm{\nabla }\cdot \left(g,\mathit{u}\right)+\mathit{\psi }-\frac{\partial }{\partial h}\left(f,g\right).$$ The terms on the right-hand side are (i) divergence of the flux of $g$, with $\mathit{u}$ being the horizontal ice current, (ii) mechanical redistribution ($\mathit{\psi })$ (i.e. ridging/rafting) and (iii) thermodynamical processes, with $f=\text{d}h/\text{d}t$ the net ice growth/melt rate. In practice, the thickness distribution is discretized over (typically 5) thickness categories (Bitz et al., 2001; Lipscomb, 2001), each characterized by a specific areal fraction (referred to as concentration). The ice thickness in each category is free to evolve between fixed boundaries.

The state of the ice is defined by a series of state variables, $X(x,y,h,t,z)$, namely ice concentration, ice volume per unit area, ice internal energy, ice salt content, snow volume per unit area and snow internal energy. Ice internal energy is the only state variable that also depends on the vertical depth in the ice ($z)$. Ice salt content does not depend on $z$ since implicit vertical salinity profiles are assumed. Following the discretization of thickness space, state variables are characterized by specific values in each category. In addition, in order to resolve the vertical profiles of internal energy, each category is further vertically divided into one layer of snow and several ice layers of equal thicknesses.

In practice, sea ice state variables follow an equation of the form $$\frac{\partial X}{\partial t}=-\mathrm{\nabla }\cdot \left(X,\mathit{u}\right)+{\mathrm{\Psi }}^X+{\mathrm{\Theta }}^X,$$ where$\mathrm{\nabla }\cdot \left(X,\mathit{u}\right)$ is the divergence of the flux of $X$, ${\mathrm{\Psi }}^X$ is the ridging/rafting and ${\mathrm{\Theta }}^X$ is the halo-thermodynamics.

Dynamics

Ice dynamics (momentum equation, advection and diffusion of state variables) is formulated on a C-grid, which is a specificity of LIM3.

Momentum equation

The ice velocity is considered the same for all categories and is determined from the two-dimensional momentum equation $$m\frac{\partial \mathit{u}}{\partial t}=A\left({\mathit{\tau }}_{\text{a}},+,{\mathit{\tau }}_{\text{w}}\right)-mf\mathit{k}\times \mathit{u}-mg\mathrm{\nabla }\mathit{\eta }+\mathrm{\nabla }\cdot \mathit{\sigma },$$ where $m$ is the ice mass per unit area, $A$ is concentration, ${\mathit{\tau }}_{\text{a}}$ and ${\mathit{\tau }}_{\text{w}}$ are the air–ice and ocean–ice stresses, $-mf\mathit{k}\times \mathit{u}$ is the Coriolis force, $-mg\mathrm{\nabla }\mathit{\eta }$ is the pressure force due to horizontal sea surface tilt and $\mathrm{\nabla }\cdot \mathit{\sigma }$ refers to internal forces arising in response to deformation. Momentum advection is at least 1 order of magnitude smaller than acceleration and is neglected (Leppäranta, 2005). The external stress terms are multiplied by concentration to satisfy free drift at low concentration (Connolley et al., 2004). The stress tensor $\mathit{\sigma }$ is computed using the C-grid EVP formulation of Bouillon et al. (2009, 2013). EVP (Hunke and Dukowicz, 1997) regularises the original VP approach (Hibler III, 1979). VP assumes a viscous ice flow (stress proportional to deformation) at very small deformations, and a plastic ice flow (stress independent of deformation) above a plastic failure threshold. This threshold lies on a so-called yield curve that depends on the ice strength determined by default from Hibler III (1979): $$P=P^{\ast }\stackrel{\mathrm{‾}}{H}{e}^{-C(\mathrm{1}-A)},$$ where $P^{*}$ and $C$ are empirical positive parameters, and $\stackrel{\mathrm{‾}}{H}$ is the ice volume per grid cell area. Other strength formulations are available in the code (e.g. Rothrock, 1975; Lipscomb et al., 2007); see Vancoppenolle et al. (2012) for details. By introducing artificial damped elastic waves and a time-dependence to the stress tensor, the EVP method enables an explicit resolution of the momentum equation with a reasonable number of sub-time steps ($\sim$ 100) and easy implementation on parallel architectures. However, EVP has to be used carefully since even the modified EVP of Bouillon et al. (2013) hardly converges to the VP solution unless a very large number ($>$ 500) of iterations is used (Kimmritz et al., 2015).

Horizontal transport and diffusion

The sea ice state variables are transported horizontally using the second-order moment-conserving scheme of Prather (1986). This scheme is weakly diffusive and preserves positivity of the transported ice fields. To smooth the ice fields and dampen instabilities, a horizontal diffusion of the form $D{\mathrm{\nabla }}^{\mathrm{2}}X$ is implemented in Eq. (3), where $D$ is a diffusion coefficient that is proportional to mean grid cell size (the reference value is 350 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ at 2${}^{\circ }$ resolution). Horizontal diffusion is solved using a Crank–Nicholson scheme, with a prescribed diffusivity within the ice pack that drops to zero at the ice edge. Horizontal diffusion should be understood as a numerical artefact introduced to avoid non-linearities arising from the coupling between ice dynamics and transport; hence, $D$ should be as small as possible.

Ridging and rafting ${\mathrm{\Psi }}^X$

The dissipation of energy associated with plastic failure under convergence and shear is accomplished by rafting (overriding of two ice plates) and ridging (breaking of an ice plate and subsequent piling of the broken ice blocks into pressure ridges). Thin ice preferentially rafts whereas thick ice preferentially ridges (Tuhkuri and Lensu, 2002). In LIM3.6, the amount of ice that rafts/ridges depends on the strain rate tensor invariants (shear and divergence) as in Flato and Hibler (1995), while the ice categories involved are determined by a participation function favouring thin ice (Lipscomb et al., 2007). The thickness of ice being deformed ($h^{\prime }$) determines whether ice rafts ($h^{\prime }<0.75$ m) or ridges ($h^{\prime }>0.75$ m), following Haapala (2000). The deformed ice thickness is $\mathrm{2}{h}^{\prime }$ after rafting, and is distributed between $\mathrm{2}{h}^{\prime }$ and $\mathrm{2}\sqrt{H^{\ast }{h}^{\prime }}$ after ridging, where $H^{*}=100$ m (Hibler III, 1980). Newly ridged ice is highly porous, effectively trapping seawater. To represent this process, mass, salt and heat are extracted from the ocean into a prescribed volume fraction (30 %) of newly ridged ice (Leppäranta et al., 1995). Hence, in contrast with other models, the net thermodynamic ice production during convergence is not zero in LIM, since mass is added to sea ice during ridging. Consequently, simulated new ridges have high temperature and salinity as observed (Høyland, 2002). A fraction of snow (50 %) falls into the ocean during deformation.

Halo-thermodynamics ${\mathrm{\Theta }}^X$

Thermodynamics refers to the processes locally affecting the ice mass and energy, and involving energy transfers through the air–ice–ocean interfaces. Halo-dynamics refers to the processes impacting sea ice salinity. In the code, both processes are assumed purely vertical and their computations are repeated for each ice category. Therefore, the reference to ice categories is implicit in this section.

Energy

The change in the vertical temperature profile, $T(z,t)$, of the snow–ice system derives from the heat diffusion equation $$\mathit{\rho }\frac{\partial E(S,T)}{\partial t}=\frac{\partial }{\partial z}\left[k,(,S,,,T,),\frac{\partial T}{\partial z}\right]+R,$$ where $z$ is the vertical (layer) coordinate, $\mathit{\rho }$ the snow/ice density (assumed constant), $E$ the snow/ice internal energy per unit mass (Schmidt et al., 2004), $S$ the salinity, $k$ the thermal conductivity (Pringle et al., 2007) and $R$ the internal solar heating rate. The effect of brine inclusions is represented through the $S$ and $T$ dependency of $E$ and $k$ (e.g. Untersteiner, 1964; Bitz and Lipscomb, 1999). The surface energy balance (flux condition) and a bottom ice temperature at the freezing point provide boundary conditions at the top and bottom interfaces, respectively. Equation (6) is non-linear and is solved iteratively. Change in ice salinity is assumed to conserve energy; hence, any salt loss implies a small temperature increase.

The solar energy is apportioned as follows. The net solar flux penetrating through the snow–ice system is $\left(\mathrm{1},-,\mathit{\alpha }\right)F_{\text{sol}}$, where $\mathit{\alpha }$ is the surface albedo and $F_{\text{sol}}$ is the incoming solar radiation flux. Only a prescribed fraction $i_{\mathrm{0}}$ of the net solar flux penetrates below the surface and attenuates exponentially, whereas the rest is absorbed by the surface where it increases the surface temperature. The radiation term in Eq. (6) derives from the absorption of the penetrating solar radiation flux $R=-\partial /\partial z\left[i_o,\left(\mathrm{1},-,\mathit{\alpha }\right),F_{\text{sw}},\exp ,\left(-,z,\mathit{\kappa }\right)\right]$, where $\mathit{\kappa }=\mathrm{1}$ m${}^{-\mathrm{1}}$ is the attenuation coefficient in sea ice, in the range of contemporary observations (Light et al., 2008). At this stage no shortwave radiation penetration is allowed when snow is present ($i_{\mathrm{0}}=\mathrm{0}$). The solar radiation flux penetrating down to the ice base is sent to the ocean. The surface albedo is a function of the ice surface temperature, ice thickness, snow depth and cloudiness (Shine and Henderson-Sellers, 1985).

Mass

The ice mass increases by (i) new ice formation in open water, (ii) congelation at the ice base, (iii) snow–ice formation at the ice surface and (iv) entrapment and freezing of seawater into newly formed ridges. It decreases by melting at both (v) the surface and (vi) the base. The snow mass increases by snowfall and reduces by surface melting, sublimation, snow–ice formation and snow loss during ridging/rafting.

Freezing and melting (i, ii, v, vi) depend on the appropriate interfacial net energy flux (open water–atmosphere, ice–atmosphere or ice–ocean) $\mathrm{\Delta }Q$ (W m${}^{-\mathrm{2}}$) such that the ocean-to-ice mass flux $F^{\text{m}}$ (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is written as $$F^{\text{m}}=\frac{\mathrm{\Delta }Q}{\mathrm{\Delta }E}.$$ $\mathrm{\Delta }E$ (J kg${}^{-\mathrm{1}}$) is the energy per unit mass required for the phase transition. For new ice formation in open water, the new ice thickness must be prescribed (usually 10 cm) and the fractional area is derived from Eq. (7). For surface melting, $\mathrm{\Delta }Q$ is different from zero only if the surface temperature is at the freezing point.

Snow-ice formation requires negative freeboard, which occurs if the snow load is large enough for the snow–ice interface to lie below sea level (Leppäranta, 1983). Seawater is assumed to flood the snow below sea level and freeze there, conserving heat and salt during the process (Fichefet and Morales Maqueda, 1997; Vancoppenolle et al., 2009b). The associated ocean-to-ice mass flux is $$F^{\text{m}}=\left({\mathit{\rho }}_{\text{i}},-,{\mathit{\rho }}_{\text{s}}\right)\frac{\partial h}{\partial t}.$$ Every ice–ocean mass exchange involves corresponding energy and salt exchanges (Schmidt et al., 2004). For instance, seawater freezing involves a change in energy $\mathrm{\Delta }E=E_{\text{i}}(S,T)-E_{\text{w}}\left(T_{\text{w}}\right)$, where $E_{\text{i}}$ is the internal energy of the frozen ice at its new temperature and salinity and $E_{\text{w}}$ is the internal energy of the source seawater at its original temperature. To ensure heat conservation in the ice–ocean system, the heat flux $Q^{\text{m}}=E_{\text{w}}\left(T_{\text{w}}\right)F^{\text{m}}$ is extracted from the ocean. Conversely, when ice melts the internal energy of melt water is sent to the ocean. Salt exchanges are detailed hereafter.

Salt

The salinity of the new ice formed in open water is determined from ice thickness, using the empirical thickness–salinity relationship of Kovacs (1996). One originality of LIM3 is that the vertically averaged ice salinity $\stackrel{\mathrm{‾}}{S}$ (in ‰) evolves in time, following Vancoppenolle et al. (2009a, b): $$\frac{\partial \stackrel{\mathrm{‾}}{S}}{\partial t}=\sum _j\left(\frac{{\mathit{\nu }}_j{S}_{\text{w}}-\stackrel{\mathrm{‾}}{S}}{h}\right)\frac{\partial h_j}{\partial t}+\sum _j{I}_j\left(\frac{S_j-\stackrel{\mathrm{‾}}{S}}{T_j}\right).$$ The first term on the right-hand side is the salt uptake summed over the three ice growth processes (ii, iii and iv), each characterized by a growth rate $\partial h_j/\partial t$ and a coefficient ${\mathit{\nu }}_j$ that determines the fraction of trapped oceanic salinity $S_{\text{w}}$. For basal freezing, ${\mathit{\nu }}_j$ is a function of growth rate (Cox and Weeks, 1988). For snow–ice formation, it is a function of snow and ice densities. For ridging, it depends on ridge porosity. The second term on the right-hand side is the salt loss summed over the two parameterized brine drainage processes (gravity drainage and flushing; see Notz and Worster, 2009). $I_j$ is 1 if the drainage process is active and 0 if it is not. Gravity drainage occurs if ice is growing at the base; flushing occurs if the snow/ice is melting at the surface. $S_j$ (5 ‰ for gravity drainage; 2 ‰ for flushing) is the restoring salinity for each drainage process and $T_j$ is the corresponding timescale (20 days for gravity drainage, 10 days for flushing).

The shape of the vertical salinity profile depends on $\stackrel{\mathrm{‾}}{S}$, so that ice with $\stackrel{\mathrm{‾}}{S}>4.5$ ‰ has a constant vertical profile. By contrast, ice fresher than this threshold has a linear profile with a lower salinity near the surface. This difference is important to properly represent the impact of brine on thermal properties (Vancoppenolle et al., 2005). Ice formation retrieves salt from the ocean, but the conjunction with water mass loss makes the ocean surface saltier. Conversely, ice melting releases salt but makes the ocean fresher. Because the ice density is assumed constant, brine drainage cannot be associated with an ice–ocean water mass exchange (the ice density would have to change to be conservative). The brine drainage flux is therefore represented as a salt flux, which directly increases ocean salinity.

Transport in thickness space

Ice growth or melt in a given category involves a transfer of ice to neighbour categories, which is formally analogous to a transport in thickness space with a velocity equal to the net growth rate $\text{d}h/\text{d}t$. This transport in thickness space is solved using the semi-Lagrangian linear remapping scheme of Lispcomb (2001). This scheme is weakly diffusive, converges with a few categories and its computational cost is minimal, which is an important property since transport operates over each ice category. Transport in thickness space is applied to all other state variables as well.

New features in LIM3.6

Control of the mass, heat and salt budgets

Mass, heat and salt must be perfectly conserved over sufficiently long timescales in an ice–ocean modelling system, especially for climate studies. Moreover, a clear identification of the different physical processes and their contributions to the air–ice–ocean exchanges is needed. These requirements were not satisfied in LIM3.0 mostly because of the temporal scheme and numerous small conservation leaks, which have necessitated a large rewriting of the code.

The changes in the sea ice state variables due to dynamics and thermodynamics were previously calculated in parallel, starting from the same initial state (Fig. 1a). Both tendencies were then combined to calculate the new state variables. This method, numerically stable and matching NEMO's philosophy, required, however, a final correction step to impose that ice losses (by melting and/or divergence) did not exceed the ice initially available. This correction step could be as important as the physical processes in some cases, and could not be attributed to a specific process. The modified temporal scheme is fractional (as for most sea ice models), removing the need for a correction step. The dynamic and thermodynamic processes are split in time and are applied sequentially (Fig. 1b), which allows for consistent diagnostics of the processes contributing to the air–ice–ocean exchanges without altering the general model behaviour (not shown). These process diagnostics are illustrated for global and regional simulations in Sects. 4 and 5.

Illustration of the changes in the time scheme. (a) The original time scheme used in LIM3.0 treats ice dynamics and thermodynamics in parallel, requiring a correction step to ensure that the ice mass is strictly positive. (b) The new scheme of version 3.6 uses an operator splitting approach, so that dynamics is calculated before thermodynamics, and therefore no correction is needed.

[Figure omitted. See PDF]

Based on these modifications, the conservation of mass, salt and heat was then carefully inspected, leading to several small corrections. In particular, the space-centred implicit backward-Euler scheme used to solve the heat diffusion equation (Eq. 6, Bitz and Lipscomb 1999) does not strictly conserve heat. The scheme is the same as in CICE, for which the problem was already reported but not yet resolved (Hunke et al., 2013). Because Eq. (6) is non-linear ($E$ and $k$ depend non-linearly on $T)$, the numerical procedure has to be iterative. The iteration stops once the temperature change is less than 10${}^{-\mathrm{5}}$ ${}^{\circ }$C or after 50 iterations. The scheme does not strictly converge, leading to an error on the heat conduction flux of $\sim$ 0.005 W m${}^{-\mathrm{2}}$, averaged over the ice pack for a global 2${}^{\circ }$-resolution simulation, with maxima reaching in some rare cases $O$ (10 W m${}^{-\mathrm{2}})$. These errors are similar to those reported in CICE user's guide (0.01 W m${}^{-\mathrm{2}}$, Hunke et al., 2013). Therefore, to ensure strict conservation, either the heat conduction fluxes or the ice temperature must be adjusted at the end of iteration. We chose to keep the ice temperature unchanged and to recalculate the net downward heat flux reaching the ocean, which could be easily implemented in other models using the same scheme.

Lateral boundary conditions

NEMO can be used in regional configurations. The BDY tool, handles the specification of boundary conditions in NEMO, with possible inflows/outflows through open boundaries (Chanut, 2005). The ocean temperature, salinity and baroclinic velocity are treated with a flow relaxation scheme (Engedahl, 1995), while the Flather (1976) radiation condition is well-suited for tidal forcing and therefore is used for both the barotropic ocean velocity and sea surface height. However, sea ice was missing from BDY, which restricted the use of regional configurations to ice-free areas. New developments to BDY were introduced to accommodate sea ice. The treatment of open boundaries in the sea ice model is not very much documented in the literature; hence, we found it difficult to compare this new approach to what is done in other models.

The sea ice state variables imposed at the boundary depend on the direction of ice velocity in a similar way to an upstream advection scheme. They are relaxed toward interior domain values where ice exits the domain, and toward external boundary data where ice enters the domain. External boundary data can either come from observations, reanalyses or reference simulations. As ice velocities in these external files are not always well determined, they need to be defined at the boundary. The tangent ice velocity is imposed to 0. The normal ice velocity depends on the presence of ice in the adjacent cell: if ice-free, ice velocity is relaxed to ocean velocity; otherwise, velocity is relaxed to the ice velocity of the adjacent cell.

Most boundary data sets do not include multiple ice categories. Hence, a strategy to fill in thickness categories in a smooth and consistent way with the external data set is defined, following the algorithm used to initialise the sea ice state variables (Vancoppenolle et al., 2012). The basic assumption relies on a distribution of ice concentration as a function of ice categories following a Gaussian law in a volume-conserving way, preserving positivity. The largest concentration is associated with the category where the mean thickness (over the grid cell) lies. Illustration of the capability of LIM3 in a regional domain is presented in Sect. 5.

Category boundaries

The original discretization of the thickness category boundaries in LIM3 follows the hyperbolic tangent formulation from CICE (Hunke et al., 2013). The formulation proved to be suitable to simulate the Arctic ice pack with only five ice categories, but cannot be easily adjusted to different ice conditions. For instance, thin ice can only be finely discretized by augmenting the number of ice categories, and de facto increasing computational cost. Multiple simulations, in particular for regional configurations, call for more flexibility without additional cpu consumption. Therefore, a new discretization was implemented that can adjust the expected mean ice thickness ($\stackrel{\mathrm{‾}}{h}$) over the domain. Category boundaries lie between 0 and $\mathrm{3}\stackrel{\mathrm{‾}}{h}$ and are determined using a fitting function proportional to $(\mathrm{1}+h{)}^{-\mathit{\alpha }}$, where $\mathit{\alpha }=0.05$. For $\stackrel{\mathrm{‾}}{h}=\mathrm{2}$ m, the new formulation is very similar to the original one. For $\stackrel{\mathrm{‾}}{h}=\mathrm{1}$ m, boundaries tighten within 3 m, providing more resolution for thin ice (Fig. 2).

Thickness category boundaries (m) as a function of categories (5 or 10). The tanh formulation from CICE, which is used in the former version 3.0 of LIM, is represented in grey and black for 5 and 10 categories, respectively. The formulation used in the new version 3.6 of LIM is proportional to $(\mathrm{1}+h{)}^{-\mathit{\alpha }}$ , where $\mathit{\alpha }=0.05$, and does not depend on the number of categories. It is displayed above for three different mean ice thicknesses $\stackrel{\mathrm{‾}}{h}$ (1, 2 and 3 m), $\stackrel{\mathrm{‾}}{h}=\mathrm{2}$ m being the closest to the tanh formulation.

[Figure omitted. See PDF]

Virtual thickness distribution

Some users may want to run LIM3.6 at the smallest possible computational cost. The most efficient way to achieve this is to use a single ice thickness category (mono-category). However, this deteriorates the results because of the poor representation of the growth and melt of thin ice, which typically reduces the amplitude of the seasonal cycle of ice extent (Holland et al., 2006; Massonnet et al., 2011). To lessen this problem, two parameterizations from LIM2 (Fichefet and Morales Maqueda, 1997) were implemented in LIM3.6. The first parameterization enhances the sea ice and snow thermal conductivities, in order to increase basal ice growth, as thin ice would do if it was properly resolved. The second parameterization aims at representing the impact of melting thin ice on ice concentration. With these two parameterizations, a mono-category simulation mostly reproduces the global mean volume and extent of a multi-category simulation, but regional differences subsist. In addition, although the mono-category approach in LIM3 is conceptually comparable to LIM2, simulations using the two sea ice models would show different results because of the different representations of halo-thermodynamics. This will be described in more details in a forthcoming contribution.

Embedded sea ice

Sea ice has been considered so far as levitating above the ocean in LIM3, and all the studies (including this one) have been based on this approximation. Even though exchanges between the levitating ice and the ocean modify the sea surface height and thermohaline structure of the ocean surface, the sea surface depression resulting from the weight of the ice and snow cover was not taken into account. The effect of embedding ice into the ocean can now be activated at will (not illustrated in this study). It improves the physical realism and influences ocean dynamics (mostly at the ice edge) via strengthened gradients of sea surface height, but does not directly affect ice dynamics (Campin et al., 2008).

A flux redistributor for the ice–atmosphere interface in coupled mode

NEMO-LIM3.6 can also be coupled to atmospheric models. Some atmospheric models only provide ice–atmosphere heat or mass fluxes ($<F>)$ for the entire grid cell, and not for each thickness category, as LIM needs. Yet the ice-atmosphere flux strongly depends on the ice surface temperature, which substantially differs among categories. To better estimate the ice–atmosphere flux in the $l$th category ($F_l)$, a “flux redistributor” has been implemented using the following linearisation: $F_l=<F>+\frac{\partial F}{\partial T^{\text{su}}}(T_l^{\text{su}}-<T^{\text{su}}>)$, where $\frac{\partial F}{\partial T^{\text{su}}}$ is the flux derivative given by the atmospheric model. $T_l^{\text{su}}$, is the ice surface temperature in the $l$th category and $<T^{\text{su}}>$ is the average over the categories, weighted by their areal fractions. The flux redistributor proves much closer to an exact computation of ice–atmosphere fluxes than a category-averaged flux.

Inputs and outputs

LIM3.6 has been interfaced with XIOS (XML input output server; http://forge.ipsl.jussieu.fr/ioserver/), a new and innovative library developed at Institut Pierre-et-Simon Laplace (IPSL) and dedicated to climate modelling data output. XIOS combines flexibility and performance. It considerably simplifies output definition and management by outsourcing output description in an external XML file. In addition, the interface offers numerous possibilities for variables manipulations such as complex temporal operations and computations involving several variables. XIOS also achieves excellent performance on massively parallel supercomputers by using several server processes exclusively dedicated to output files. File system writing is performed concurrently with computation.

Global ice–ocean simulation: ORCA2-LIM3

Experimental set-up and observation data sets

The simulation presented here is the standard simulation that can be performed with the most recent 3.6 version of NEMO right after downloading the code, in one of the main supported NEMO configurations (ORCA2-LIM3), and forced by the reference CORE normal year forcing directly provided with the code. This is not the best simulation that can be produced, but rather the one that a user starting with the model would perform.

In ORCA2-LIM3, NEMO comprises the ocean general circulation model OPA version 3.6 (Madec, 2008) and LIM (Vancoppenolle et al., 2009a) in its 3.6 version presented above, running on the same 2${}^{\circ }$-resolution grid (ORCA2). More details can be found in Mignot et al. (2013). The atmospheric state is imposed using the CORE normal year forcing set proposed by Large and Yeager (2009), developed to inter-compare ice–ocean models (e.g. Griffies et al., 2009). It is based on a combination of NCEP/NCAR reanalyses (for wind, temperature and humidity) and various satellite products (for radiation), has a 2${}^{\circ }$ resolution and near-zero global mean heat and freshwater fluxes. The so-called normal year data set superimposes the 1995 synoptic variability on the mean 1984–2000 seasonal cycle. The simulation lasts 100 years, much longer than needed for sea ice to reach equilibrium. Most diagnostics presented hereafter are seasonal averages over the last 10 years of the simulation. The computational cost of such a simulation is about 12 h on 64 processors of an IBM Power6, with LIM3.6 consuming less than 25 % of this time.

Mean sea ice concentrations from the simulation ORCA2-LIM3 and the observations OSI-SAF for March and September in the Arctic (left panels) and February and September in the Antarctic (right panels). The white line indicates the 15 % ice concentration contour.

[Figure omitted. See PDF]

The observed ice extent is derived from ice concentration retrievals of the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI-SAF; Eastwood et al., 2010) and is presented here as 1984–2000 monthly means. To put the simulated ice volume in context, we do not use satellite estimates, for which uncertainties are very large (e.g. Zygmuntowska et al., 2014), but instead the 1979–2011 reanalysis PIOMAS in the Arctic (Schweiger et al., 2011), and the NEMO-LIM2-EnKF reconstruction in the Antarctic (Massonnet et al., 2013).

Ice concentration and thickness

Neither the model nor the atmospheric forcing are precisely tuned to get the most realistic sea ice simulation, because this depends on forcing, resolution and user wishes. Instead, we choose the model default parameters with the standard reference forcing and show that the simulated ice concentrations and thicknesses are in reasonable agreement with observations.

Mean seasonal cycle of sea ice extent (i.e. area inside the 15 % concentration contour) in the Northern (in blue) and Southern (in cyan) hemispheres from the ORCA2-LIM3 simulation (solid lines) and as derived from OSI-SAF observations (dashed lines). Units are in 10${}^{\mathrm{6}}$ km${}^{\mathrm{2}}$.

[Figure omitted. See PDF]

Figure 3 shows the ice concentrations at the model maximum and minimum extent in ORCA2-LIM3 and OSI-SAF (March and September for the Arctic; February and September for the Antarctic). The simulated ice distribution is relatively close to the observations, with some common defects. In the boreal winter, the ice extends too much southward covering a large part of the Greenland Sea, while it almost disappears near Antarctica. These biases have unclear origins and we do not intend to resolve them but some leads can be proposed. In the Northern Hemisphere, we notice a low ocean heat supply by the North Atlantic Current and an overestimated ice volume export through Fram Strait, which could explain some of the bias. But other factors as the forcing or model physics, in particular dynamics, cannot be ruled out. In the Southern Hemisphere, we notice a wrong position of the Antarctic Circumpolar Current and an overestimated ocean convective activity, which melts ice by mixing relatively warm and salty water at depth with cold and fresh surface waters, and which could explain the ice loss. Such problems are common in global ocean models (Kim and Stössel, 2001), and vertical physics in the ocean should certainly be tuned to improve the realism of the simulated ice characteristics.

Mean simulated sea ice thicknesses (m) at the time of maximum ice volume: for March in the Northern Hemisphere and for September in the Southern Hemisphere.

[Figure omitted. See PDF]

The seasonal cycle of the sea ice extent (i.e. the area enclosed within the 15 % ice concentration contour, white lines in Fig. 3) is presented in Fig. 4 for both hemispheres. The model reproduces the amplitude of the observed seasonal variations of ice extent but is biased low all year long, and especially in austral summer.

The simulated ice thickness distributions are displayed in Fig. 5 for both hemispheres, at the time of maximum extent (March and September). The ice thickness exceeds 3 m in the central Arctic, reaching 5 m along the Canadian and Greenland coasts. This is in rough agreement with the submarine thickness retrievals (3.4 m in the central Arctic in February–March 1988; Kwok and Rothrock, 2009). The spatial distribution follows expectations, except a spurious band of thick ice along the East Siberian shelf. The simulated Arctic ice volume ranges from 17 000 km${}^{\mathrm{3}}$ in September to 35 000 km${}^{\mathrm{3}}$ in March–April, i.e. somewhat higher than Pan-arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) reanalyses. In the Southern Hemisphere, the ice is generally thinner than in the Arctic, with a modal value of nearly 1 m. The model underestimates the thickness of thick ice in the Weddell and Amundsen seas (Worby et al., 2008; Kurtz and Markus, 2012). The band of thick ice along the east side of the Antarctic Peninsula is missing, which is attributed to misrepresented NCEP winds in the region (Timmermann et al., 2005; Vancoppenolle et al., 2009b). The simulated ice volume (0 to 14 000 km${}^{\mathrm{3}})$ is somewhat larger than the reanalysis values (2000–10 000 km${}^{\mathrm{3}}$; F. Massonnet personal communication, 2015) and satellite estimates (3000–11 000 km${}^{\mathrm{3}}$, Kurtz and Markus, 2012).

Simulated mean seasonal cycles of the different ice mass balance processes in the ORCA2-LIM3 simulation: Arctic (left panel) and Antarctic (right panel). Ice grows from the base (magenta), in open water (red), by snow–ice formation (orange) or by freezing of sea water trapped in the ridges (green). Ice melts at the base (blue) and surface (cyan). Ice advection is nil here since diagnostics are hemispheric. The black line is the net ice production (i.e. the sum of all the processes). Units are in cm month${}^{-\mathrm{1}}$. Positive and negative values represent creation and destruction of sea ice, respectively.

[Figure omitted. See PDF]

This simulation could obviously be improved through careful calibration, which depends on resolution and forcing. Calibration can be achieved by adjusting the atmospheric forcing and vertical ocean physics, and by tuning the most influential ice parameters. For instance, the Arctic ice thickness can be increased substantially by increasing the albedo, decreasing the minimum lead fraction or decreasing ice strength.

Mass and salt balances

The new developments allow for an examination of the ice mass, heat and salt budgets seasonally and over the different processes. Seven processes affect the ice mass (see Sect. 2.3.2). Five belong to vertical thermodynamics: new ice growth in open water, basal growth and melt, surface melt and snow–ice formation. Two are dynamical processes: advection and entrapment and freezing of seawater in newly built ridges. Changes in the heat and salt contents involve the same processes, plus the changes in internal temperature (for heat budget) and internal salinity due to brine drainage (for salt budget).

We focus on the mass budget for illustration and present its different contributors integrated over the Northern and Southern hemispheres in Fig. 6. In both hemispheres, the dominant balance is between basal ice growth and melt. Surface melting is also important but only in the Arctic during boreal summer. Contributions of secondary importance are new ice formation in open water during the cold season (both hemispheres) and snow–ice formation during Antarctic spring. Note that the contribution from advection is obviously nil when integrated over a hemisphere. The maximum growth rate is about the same in both hemispheres (slightly larger than 20 cm month${}^{-\mathrm{1}}$). Basal melt is remarkably weaker in the Arctic than in the Antarctic (maximum at 40 and 70 cm month${}^{-\mathrm{1}}$, respectively). This is because in the Arctic, the ice is constrained by continents to stay at high latitudes, where the ocean stratification is strong and the ocean heat flux is weak. Overall, about 26 000 km${}^{\mathrm{3}}$ of ice are formed and melted each year in the Arctic, which corresponds to about 2 m of ice. About 320 Gt of salt are extracted from the ocean during freezing and released during ice desalinisation and melting. These mean values are similar in the Antarctic: 22 000 km${}^{\mathrm{3}}$ of annual ice production ($\sim$ 1.8 m) and 320 Gt of salt.

This integrated view masks strong geographic disparities. In Fig. 7 we show the geographical distribution of some of the processes in March in the Arctic. The interior of the ice pack still grows from the bottom, while the base of the ice edge melts, resulting in snow–ice formation where snow is thick enough. As expected, the strongest thickness changes due to advection are near the ice edge. Ice formation in open water is globally weak but becomes one of the main processes in some regions of climate importance (see next section).

Horizontal distribution of the five relevant processes contributing to the sea ice mass balance in March in the Northern Hemisphere, from the ORCA2-LIM3 simulation. Units are in cm day${}^{-\mathrm{1}}$. Positive and negative values represent creation and destruction of sea ice, respectively.

[Figure omitted. See PDF]

Regional configurations

Experimental set-up

To illustrate the capability of NEMO-LIM3 in regional ice-covered domains, we designed an experiment in a regional configuration (500 $\times$ 500 km) around the Svalbard Archipelago. This region was chosen because of the diverse conditions encountered and strong tides (a tidal gauge at Ny-Ålesund, on the west coast of Svalbard, records tidal amplitudes up to 2 m). North of the archipelago, lies the perennial ice pack of the Arctic Ocean transitioning southwards to a seasonal ice zone. The domain also includes the large Storfjorden polynya, frequently open during winter. Polynyas are small (10–10${}^{\mathrm{5}}$ km${}^{\mathrm{2}})$ and sporadic by nature, but their role in climate is important (e.g. Morales Maqueda et al., 2004). In winter, the ocean heat loss in polynyas is considerable, producing large amounts of sea ice, as well as dense water sinking towards the deep ocean basins. At the onset of melting season, polynyas enhance ice melting as the open waters capture more heat than ice-covered areas.

Horizontal resolution is very high (2 km) in order to properly represent fine-scale processes taking place in polynyas. The basin is vertically discretized by 75 non-uniform ocean levels, with a resolution of 1 m at the surface. The domain is open at the four boundaries and conditions there are set up using the BDY tool, modified as described in Sect. 3.2. Bathymetry is interpolated from etopo1 (Amante and Eakins, 2009), which actually retrieves data from IBCAOv2 north of 64${}^{\circ }$ N (Jakobsson et al., 2008). Tides are important drivers for high-frequency processes. Therefore, they are included here as well as the non-linear free surface ($z^{*}$ coordinates system). A third-order upstream biased advection scheme is used for ocean tracers and momentum (instead of the flux corrected transport used in ORCA2-LIM3). Such a scheme is indeed more precise and has implicit diffusion. It also minimizes diffusion; hence, the oceanic structures can develop without being impeded by homogeneous diffusion. The $k$-$\mathit{\epsilon }$ closure scheme using generic length-scale turbulent mixing is chosen (Umlauf and Burchard, 2003; Reffray et al., 2015). The simulation is forced at the surface by a 6-hourly, 3/4${}^{\circ }$ $\times$ 3/4${}^{\circ }$ ERAI data set, and at the boundaries by 5-day outputs from a DRAKKAR 1/4${}^{\circ }$ global reference simulation ORCA025-MJM (an update to ORCA025-G70; Barnier et al., 2006; Molines et al., 2007; DRAKKAR group, 2007). We also prescribe tidal sea surface height and barotropic velocity at the boundaries from FES2012 (Carrère et al., 2012). The simulation is conducted over 1999–2009 in order to capture inter-annual variability.

The model behaviour at the boundary is satisfactory. No noise or wave reflection pollutes the basin despite strong in and out flows and the presence of tides (not shown). The simulation is also able to represent transient polynya occurrence. As an example, Fig. 8 shows the simulated ice concentrations on 22 May 2002 around Svalbard (right panel) as well as the corresponding observations (left panel). At this date, northeastern winds were sufficiently strong to open the Storfjorden polynya by pushing sea ice towards the western side of the fjord. The simulated opening of polynyas – in terms of timing, location and size – is reasonable in Storfjorden and elsewhere, though polynyas are somewhat smaller than observed and their location is not precisely captured. This is likely due to the low spatio-temporal variability of the ERAI surface forcing, as highlighted by previous studies (Skogseth et al., 2007). Downscaling the forcing with a regional atmospheric model is probably required to further improve the simulation. The Storfjorden polynya is not exactly found where it should be, north of Storfjorden (Fig. 8), which could be due to the atmospheric forcing or to the absence of a representation of landfast ice in the model and must be further investigated.

(Left)satellite MODIS image of the Svalbard Archipelago (22 May 2002). Note that clouds and sea ice are both white. (Right) 1-day averaged simulated sea ice concentrations at the same date from the high-resolution regional simulation. In both pictures, ice is pushed away from the shore by northeasterly winds, allowing formation of a polynya along the east coast of Storfjorden.

[Figure omitted. See PDF]

Mass and salt balances in Storfjorden

Figure 9 shows the 10-year variability of the different mass balance processes over the 13 000 km${}^{\mathrm{2}}$ of the Storfjorden region (see Fig. 8). The sea ice mass balance is dominated by basal growth (16 km${}^{\mathrm{3}}$ year${}^{-\mathrm{1}}$) and new ice growth in open water (12 km${}^{\mathrm{3}}$ year${}^{-\mathrm{1}}$), compensated by export out of the domain (not shown) and basal melt (11 km${}^{\mathrm{3}}$ year${}^{-\mathrm{1}}$). Surface melt can be significant (up to one third of total melt) but only at the beginning of summer. As expected, ice growth in open water is a crucial process here, while it is weak once averaged over the Arctic basin (see previous section). The net ice production is $+$17 km${}^{\mathrm{3}}$ year${}^{-\mathrm{1}}$ on average, with strong inter-annual variability (from 23 km${}^{\mathrm{3}}$ in 2001–2002 to 10 km${}^{\mathrm{3}}$ in 2006–2007). This corresponds to a salt flux from the ocean to the ice of about 150 Mt year${}^{-\mathrm{1}}$. Over a year, the net production almost balances ice export (not shown), so there is no long-term accumulation of ice in the basin. However, at timescales shorter than a year, ice can pile up in the Storfjorden.

A 10-year inter-annual variability of the processes involved in ice evolution integrated over the Storfjorden area from the regional simulation. Processes are the same as in Fig. 6, plus an advection term corresponding to ice coming in and out of the area, which is not shown for more clarity. Units are in cm day${}^{-\mathrm{1}}$. Positive and negative values represent creation and destruction of sea ice, respectively. Note that for more readability, variations are smoothed with a Hanning filter at a period of 2 months.

[Figure omitted. See PDF]

By combining AMSR-E sea ice concentrations and atmospheric forcing from ERA-interim, Jardon et al. (2014) estimated a mean ice production of 47 km${}^{\mathrm{3}}$ in winter between 2002 and 2011. With a similar approach, Skogseth et al. (2004) found a mean ice production of 40 km${}^{\mathrm{3}}$ during 1998–2002. In our simulation, this production amounts to 33 km${}^{\mathrm{3}}$ for the period 1999–2009. This value is reasonable, though it is smaller than observational retrievals and reanalyses. This could be related to the small size of the simulated polynya and/or to the lack of high-resolution, high-frequency winds in the ERAI forcing and should be further investigated.

Conclusions

The Louvain-la-Neuve sea ice model (LIM) has evolved considerably during the past decade. Two versions have been developed and have coexisted up until now. LIM2 is based on a Hibler III (1979) mono-category approach, and was integrated in the NEMO system about 1 decade ago (Timmermann et al., 2005). It was the reference model to date and was used in a variety of simulations including CMIP5. LIM3 is a more sophisticated model developed 5 years ago (Vancoppenolle et al., 2009a), including a better representation of sub-grid-scale ice thickness distribution and salinity processes. Several modifications to LIM3 have been done recently to make it more robust, versatile and sophisticated, leading to LIM3.6, described in this paper. LIM3.6 is the reference model for the forthcoming CMIP6 simulations, while LIM2 is no longer the reference and will be discontinued in the next NEMO release.

LIM3 has been improved for a use in various configurations, from climate to regional studies, with a large range of resolutions and complexities. Three main developments were required. First, the code has been made strictly conservative. For that purpose, the general time stepping has changed from parallel to a splitting approach. In other words, thermodynamics processes are now performed after dynamics, which enables the discrimination of the different processes contributing to the mass, heat and salt exchanges across the interfaces between air, ice and ocean. Conservation in the code has been carefully examined by comparing these exchanges with thermodynamical and dynamical ice evolution, which has led to several small corrections to reach a strictly conservative code. In particular, the iterative procedure to solve the heat diffusion equation (Eq. 6) did not exactly converge, leading to small heat leaks. The leaks are now corrected by recalculating heat fluxes. Second, version 3.6 of LIM is the first to handle open-boundary conditions for regional simulations in ice-covered areas. The sea ice state variables at the boundary depends on the direction of the normal ice velocity to allow realistic inflows and outflows with the rest of the ocean. Boundary conditions are flexible enough so that ice boundary data sets can either integrate a sub-grid-scale ice thickness distribution or not. In addition, the formulation of the discretization of ice categories boundaries has changed to adapt a simulation to different ice thickness conditions, as encountered in regional configurations. Third, LIM3.6 sophistication and versatility have further increased. A mono-category capability has been implemented with the parameterization of thin ice melting, especially for users needing an ice model at minimal computational cost. A flux redistributor at the top of the ice categories has been coded for the coupling with atmospheric models that cannot handle multiple fluxes over a grid cell. Finally, the effect of the ice and snow weight on the sea surface height has been implemented.

To illustrate some of the new capabilities of LIM3, we present 100 years of the 2${}^{\circ }$-resolution forced simulation ORCA2-LIM3, and 10 years of a regional simulation at 2 km resolution around the Svalbard Archipelago, which hosts the recurrent Storfjorden polynya. We mainly focus on the ice mass budget and show how they differ, depending on the region studied. At the global scale, the dominant processes are basal ice growth and basal ice melt for both hemispheres, but other processes matter locally. In the Storfjorden, new ice growth in open water is nearly as large as basal growth. The entire ice production is exported out of Storfjorden annually. Production presents large inter-annual variability over the 10 years of the experiment (1999–2009), with maximum values exceeding twice the minimum.

There are also ongoing and upcoming developments for LIM.

• The compatibility between the Adaptive Grid Refinement In Fortran (AGRIF; Debreu et al., 2008) and LIM3 to run global simulations is yet to be achieved and work is in progress to use LIM2–AGRIF interface (Talandier et al., 2014) and apply it to LIM3.

• The melt pond parameterization of Flocco and Feltham (2007), as implemented by Lecomte et al. (2015), exists in a branch of the code and is expected soon in the reference version.

In the future, LIM will continue to be developed, including, among others, sea ice biogeochemistry (Vancoppenolle and Tedesco, 2015; Moreau et al., 2015), an elasto-brittle rheology (Girard et al., 2011), improved snow physics (Lecomte et al., 2013, 2015) and a sub-grid-scale representation of ice–ocean exchanges (Barthélemy et al., 2015).

Code availability

The version 3.6 of LIM3 is incorporated in the reference version of NEMO (currently v3.6 stable) and can be downloaded from the NEMO web site (http://www.nemo-ocean.eu/) at this address: https://forge.ipsl.jussieu.fr/nemo/browser/branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3.

Acknowledgements

We thank the two anonymous reviewers for their thorough reading and constructive comments, which helped improving the manuscript. The present work was supported by the French National Research Agency (ANR) as part of the OPTIMISM project (ANR-09-BLAN-0227-01), the European Community's Seventh Framework Program FP7/SPACE (MyOcean2, grant 283367), IS-ENES2 (grant 312979) and BISICLO (FP7 CIG 321938). Computational resources have been provided by the French Institut du Développement et des Ressources en Informatique Scientifique (IDRIS). The Antarctic sea ice concentration reprocessed data set was provided by the Ocean and Sea Ice Satellite Application Facilities (OSISAF) and is available at http://osisaf.met.no. The MODIS image is courtesy of the online Data Pool at the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://lpdaac.usgs.gov/data_access). Edited by: P. Huybrechts