2.1 Ocean–sea ice model (NEMO3.6-LIM3)

This study is conducted using the dynamic–thermodynamic sea ice model LIM3.6 (Louvain-la-Neuve sea Ice Model; Rousset et al., 2015). It is coupled to a finite-difference, hydrostatic, free-surface, primitive-equation ocean model within version 3.6 of the NEMO framework (Nucleus for European Modelling of the Ocean; Madec, 2008). LIM includes an ITD, which is described in more details in the next section, energy conserving thermodynamics (Bitz and Lipscomb, 1999), a time-dependent salinity profile (Vancoppenolle et al., 2009a) and an elastic–viscous–plastic rheology formulated on a C grid to model the ice dynamics (Bouillon et al., 2013).

2.2 Ice thickness distribution

One of the key features of LIM is the ITD that is used to represent the subgrid-scale distribution of sea ice thickness, enthalpy and salinity (Thorndike et al., 1975). The probability density function for thickness h, $g(h,\mathit{x},t)$, evolves according to the following equation:

The first term of the right-hand side is the convergence of the horizontal flux of g with u the ice velocity. In the second term, f is the thermodynamic growth rate. The third term, ψ, is the mechanical redistribution. The way this term is modeled is somewhat specific to LIM: it includes ridging and lead opening terms (as in other models) but also a rafting term (Vancoppenolle et al., 2009b). Another specificity of LIM is that newly formed ridges have a prescribed 30 % porosity.

The numerical formulation of the ITD in the model follows Lipscomb (2001). The thickness distribution is discretized into a fixed number of categories. Sea ice in each category occupies a varying areal fraction of the grid cell. Ice thickness in each category is constrained to remain between user-prescribed boundaries. Ice growth and melt induce transfers of ice volume and area between categories, which is dealt with the linear remapping scheme of Lipscomb (2001). This scheme combines the advantages of being computationally inexpensive and only weakly diffusive.

In the current version of LIM, the ITD is discretized by default by using a fitting function that places the category boundaries between 0 and $\mathrm{3}\stackrel{\mathrm{‾}}{h}$, with $\stackrel{\mathrm{‾}}{h}$ the expected mean ice thickness over the domain. A greater resolution is placed for thin ice (Rousset et al., 2015). More specifically, the lower and upper limits for ice thickness in category $i=\mathrm{1},\mathrm{\dots },N$ are f(i−1) and f(i) with

where N is the number of categories, $i=\mathrm{1},\mathrm{\dots },N$ refers to the category index and α=0.05. The upper limit of the last category is always reset from $\mathrm{3}\stackrel{\mathrm{‾}}{h}$ to 99.0 m, in order to allow hosting very thick ice.

2.3 Model configuration

The configuration of NEMO3.6-LIM3 is almost identical to the one used in Barthélemy et al. (2017), where the reader may find the details that are not recalled here. The revision 6631 of the branch 2015/nemo_v3_6_STABLE of the NEMO SVN repository (https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-3.6/, last access: 19 August 2019) is used. The ocean and sea ice models are run on the global eORCA1 grid, which has a nominal resolution of 1^∘ in the zonal direction. In the ocean, a partial step z coordinate with 75 levels is used along the vertical. The layer thickness increases non-uniformly from 1 m at the surface to 10 m at 100 m depth and reaches 200 m at the bottom. To avoid spurious model drift, a weak restoring towards the World Ocean Atlas 2013 surface salinity climatology (Zweng et al., 2013) is applied, with strength 167 mm d⁻¹ (i.e., a relaxation timescale of 10 months for a 50 m deep oceanic mixed layer). In order to avoid adversely altering ocean–ice interactions, this restoring is damped under sea ice (multiplied by 1 minus sea ice concentration), i.e., where the observations are less reliable.

Figure 1Ice thickness category boundaries in the three sets of sensitivity experiments. The upper boundary of the last category is always set to 99.0 m. Note that the ice thickness scale is different in the three panels. Because the ITD discretization in the third set of experiments (S3) branches from experiment S2.09 of the second set, that experiment is repeated in the list but labeled as S3.09.

Download

Figure 2Mean seasonal cycles of Arctic (a, b) and Antarctic (c, d) sea ice extents (a, c) and volumes (b, d), over 1995–2014, in the first set of sensitivity experiments. Ice extents derived from the Ocean and Sea Ice Satellite Application Facility (OSISAF) sea ice concentration observational product (OSI-409a; EUMETSAT, 2015) are also shown, as well as Arctic and Antarctic ice volumes derived from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) and Global Ice Ocean Modeling and Assimilation System (GIOMAS) reanalyses, respectively (Schweiger et al., 2011; Zhang and Rothrock, 2003). The stars show the monthly data and the curves are cubic interpolations between the data points.

Download

The atmospheric forcing is provided by the DRAKKAR Forcing Set version 5.2 (DFS5.2; Brodeau et al., 2010; Dussin et al., 2016). This global forcing set is derived from the ERA-Interim reanalysis over 1979–2015 and from ERA-40 for the period 1958–1978. It has a spatial resolution close to 0.7^∘, or 80 km, and is utilized within the CORE forcing formulation of NEMO, which uses bulk formulas developed by Large and Yeager (2004). Continental freshwater inputs consist of river runoff rates from the climatological dataset of Dai and Trenberth (2002) north of 60^∘ S, of prescribed meltwater fluxes from ice shelves along the coastline of Antarctica (Depoorter et al., 2013) and of climatological freshwater fluxes from iceberg melt at the surface of the Southern Ocean (Merino et al., 2016). The sole difference compared to Barthélemy et al. (2017)'s setup is the shorter integration time (1979–2015 instead of 1958–2015) due to the computational cost of the multiple sensitivity experiments. Simulations are only weakly impacted by the length of integrations. A comparison of the reference simulation labeled S1.05 (see below) with the corresponding experiment in the aforementioned study indicates that the interannual variability of hemispheric ice volumes is the only sea ice index showing a weak sensitivity to the start date, while other ones are not affected (not shown). The ocean is initialized at rest with temperature and salinity fields from the World Ocean Atlas 2013 (Locarnini et al., 2013; Zweng et al., 2013) and we analyze the last 20 years of the simulations (1995–2014).

Figure 3Effective sea ice thickness (grid-cell ice volume divided by grid-cell area) in the Arctic in March (top) and in the Antarctic in September (bottom), averaged over 1995–2014. Experiment S1.05 and differences with selected experiments of the first set are shown.

2.4 Sensitivity experiments

An ITD discretization features two characteristics: the number of categories used and the position of category boundaries. By default, the LIM automatically sets the position of category boundaries for a given number of categories, according to Eq. (2). However, this choice can be overridden by the user, in order to explore specific discretizations.

We ran three successive sets of sensitivity experiments. In the first set, labeled S1, the default ITD discretization of LIM is used and both the number and positions of category boundaries change according to Eq. (2) for a number of 1, 3, 5, 10, 30 and 50 categories. The value for $\stackrel{\mathrm{‾}}{h}$ is set to 2.0 m and represents a tradeoff between gross estimates of Antarctic and Arctic basin-wide average thickness (∼1.0 and 3.0 m, respectively). The category boundaries are displayed in the top panel of Fig. 1.

The results derived from this first set of experiments can potentially be useful for users of the model but not necessarily for understanding the physical processes driving the model response. Indeed, in S1, both the position and the resolution of the thickness categories vary from one experiment to the next; it is thus impossible to disentangle one effect from the other. Therefore, in a second set of experiments, labeled S2, we bypass the default formulation of LIM and successively append new thickness categories without changing the existing category boundaries (Fig. 1, middle panel). This set of experiments allows testing the influence of thick ice categories, as Hunke (2014) suggested this aspect to be potentially important. In that set of experiments, the resolution of the ITD is unchanged and we can test the specific role of the category positions on the simulated mean sea ice state.

Finally, in order to test the influence of the number of ice categories, independently of their position, we run a third set of sensitivity experiments (S3), in which the lower boundary of the thickest category is locked to 4.0 m as in experiment S2.09, and the ITD is coarsened/refined by merging/splitting the existing categories (Fig. 1, bottom panel). The upper limit of 4.0 m corresponds to the maximum thickness that thermodynamic ice growth can sustain in the Arctic (Maykut and Untersteiner, 1971) and allows therefore the thickest category to host the deformed ice produced by the model.

3.1 Influence of the number of categories in the default formulation (S1 experiments)

Seasonal cycles of Arctic and Antarctic sea ice extents and volumes are presented in Fig. 2. Following the usual conventions, sea ice extent is defined as the sum of all oceanic grid-cell areas that contain at least 15 % of ice concentration, while sea ice volume is defined as the sum of individual grid-cell volumes. The main consequence of using a larger number of ice thickness categories is to increase the ice volume in winter. In the Arctic, the increase persists during the whole seasonal cycle, even though it becomes smaller in summer. To position the volume produced by our model, a comparison is made with the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) reanalysis (Schweiger et al., 2011). The model produces higher volume than the PIOMAS reanalysis for simulations with more than three categories. In the Antarctic, the increase is limited to the ice growing season, while the rest of the year features a decrease in volume when using only few categories (S1.01 and S1.03), which is due to an excessive sea ice extent in summer (discussed later). In either hemisphere, the annual maximum of ice volume does not converge to an asymptotic value when increasing the number of categories. The possible origins of this lack of convergence will be discussed in the next section.

Figure 4Contributions to the seasonal mass balance of Arctic (a) and Antarctic (b) sea ice as simulated in the first set of sensitivity experiments (S1) (varying number of categories with the default LIM3 ITD discretization), averaged over the areas depicted in Fig. A1. Blue colors refer to processes that contribute to positive ice volume changes, while red colors indicate processes implying negative ice volume changes. The name of the experiments is indicated in the upper panel for the January–February–March season.

Download

3.1.1 Response of sea ice thickness

The impact of changing the number of categories on the spatial distribution of ice thickness is shown in Fig. 3. Compared to experiment S1.05, using only one category reduces the ice thickness mainly in the thick ice areas, by up to 1 m north of the Canadian Arctic Archipelago and the Canadian basin in the Arctic, and by around 0.2 m in the Ross and Weddell seas in the Antarctic. By contrast, increasing the number of categories leads to a more uniform thickening of the ice, by up to 0.5 m (0.2 m) for S1.50 in the Arctic (Antarctic).

The origin of a thickness increase with the number of ice thickness categories can be explored with the tendency diagnostics provided in LIM. Indeed, variations in state variables, including volume in each category and therefore the aggregate volume, can be attributed to various physical processes accounted for in the model such as open-water ice formation, bottom growth, snow-ice production and dynamic production (i.e., ice formed due to refreezing of seawater after entrapment into porous ridges) but also bottom melt and surface melt.

The increase in sea ice thickness is mainly due to an enhancement of thermodynamic basal growth rates in fall and winter (Fig. 4). Our hypothesis is that, for the same grid-cell average thickness, a better-resolved ITD results in larger basal ice growth rates due to the inverse relationship between conductive heat fluxes and sea ice thickness. To illustrate this point, let us consider two counterfactual configurations A and B, with the same grid-cell average ice thickness but different ITDs. In A, a uniform 2 m thick slab of ice covers 100 % of the grid cell, while in B 50 % of the grid cell is covered by a 1 m thick slab of ice and the other 50 % by a 3 m thick slab of ice. All other things being equal (surface and bottom ice temperatures, snow, ice conductivity, ocean–ice heat flux), the heat conduction flux will be proportional to 1∕2 W m⁻² in A as per Fourier's law of conduction. By contrast, the heat conduction flux will be proportional to 1 and 1∕3 W m⁻² on each half of the grid cell in B, giving an average of 2∕3 W m⁻². The basal ice growth, resulting from the imbalance between ocean–ice heat flux and the heat conduction flux, will therefore be higher in B.

Figure 5Sea ice concentration in the Arctic in August (top) and in the Antarctic in February (bottom), averaged over 1995–2014. Experiment S1.05 and differences with selected experiments of the first set are shown.

The diagnosed sea ice mass balance analysis of Fig. 4 also shows that, in the Arctic, the run with one thickness category compensates the deficit of basal ice growth (relative to multi-category runs) by enhanced dynamic production during fall and winter (October to March). We can understand this finding as follows. First, the deficit of ice growth leads to a thinner ice pack. However, in the model, compressive ice strength is a function of the grid-cell average thickness H and the concentration A:

$$\begin{array}{}\text{(3)}& P=P^{*}H{e}^{-C(\mathrm{1}-A)},\end{array}$$

where $P^{*}=\mathrm{20}$ kN m⁻² and C=20 are two empirical constants. In all simulations, A∼1 in winter due to the thermal constraint imposed by the atmospheric forcing. Small values of H in the one-category simulation lead to less resistance of ice to compression and hence enhanced mechanical redistribution, which fosters dynamic growth.

What is the origin of non-convergence of sea ice volumes noticed in Fig. 2? Based on the mass balance analysis, the non-convergence is primarily due to a non-convergence in basal growth rates (Fig. 4). However, in an idealized configuration with no snow and without ice dynamics, growth rates reach 95 % of the asymptotic value when five categories or more are used (Fig. 6). Thus, the non-convergence has to be the consequence of an interplay between the dynamic and thermodynamic processes: ridging or rafting transfers thin ice to thicker categories, allowing more ice production by thermodynamic growth. It is not possible to produce deeper analyses at this stage as the ice mass balance terms are only available at the grid-cell level, not at the ice thickness category level. It is also unclear if the non-convergence is the consequence of our experimental setup in which the atmospheric forcing is prescribed – offering no possibility for negative feedbacks to operate or at least not as strongly as they might do in a coupled model. We will soon be testing this hypothesis with a coupled model.

Figure 6(a) A supposed ice thickness distribution in a model grid cell (red; log-normal distribution with mean 3 m and standard deviation 2 m) and its discretization in five categories following the default formulation of LIM3 (grey bars). (b) Average basal growth rates for $\mathrm{1},\mathrm{2},\mathrm{\dots }\mathrm{30}$ categories. For each category, the basal growth rate was computed assuming sea ice thickness equal to the category mean (h_i with $i=\mathrm{1},\mathrm{\dots },N$, where N is the number of categories), assuming no snow, an atmosphere–ocean temperature difference of ΔT=30 K, sea ice conductivity of k=2 W mK⁻¹, latent heat of fusion of L_f=334 000 J kg⁻¹ and sea ice density of ρ_i=917 kg m⁻³. The growth rates were then averaged over categories, taking into account the relative area of each category: $\dot{h}=\frac{\mathrm{1}}{{\mathit{\rho }}_i{L}_{\mathrm{f}}}{\sum }_{i=\mathrm{1}}^N\frac{k\mathrm{\Delta }T}{h_i}g\left(h_i\right)\mathrm{d}{h}_i$, with dh_i the width of category i.

Download

3.2 Influence of thick ice categories (S2 experiments)

The aim of the second set of sensitivity experiments is to examine specifically the importance of thick ice categories in the ITD, without refining the discretization in the thin range (Fig. 1, middle panel). In this set, the sea ice volume is again the variable that is most impacted by the selected number of categories (Fig. 7). The Arctic ice volume increases monotonically with the addition of new thick categories until a threshold at ∼4 m is reached: adding categories beyond that threshold (experiments from S2.09 on) does not have an impact on the total volume. The Antarctic ice volume is less sensitive.

Figure 7Mean seasonal cycles of Arctic (a) and Antarctic (b) sea ice volumes, over 1995–2014, in the second set of sensitivity experiments. Also shown are the Arctic and Antarctic ice volumes derived from the PIOMAS and GIOMAS reanalyses, respectively (Schweiger et al., 2011; Zhang and Rothrock, 2003). The stars show the monthly data and the curves are cubic interpolations between the data points.

Download

Figure 8Mean ice thickness distribution in the second set of sensitivity experiments. For each ice thickness category, the relative areal proportion of ice for that category was estimated from the Arctic (March, top) and Antarctic (September, bottom) 1995–2014 average sea ice concentration field restricted to the domain shown in Fig. A1. Thin vertical lines delimit the category boundaries. Note that, for the sake of readability, the spacing along the x axis is logarithmic. The upper bound of the last category is always set to 99.0 m and is not displayed.

Download

Figure 9Mean seasonal cycles of Arctic (a) and Antarctic (b) sea ice volumes, over 1995–2014, in the third set of sensitivity experiments. Also shown are the Arctic and Antarctic ice volumes derived from the PIOMAS and GIOMAS reanalyses, respectively (Schweiger et al., 2011; Zhang and Rothrock, 2003). The stars show the monthly data and the curves are cubic interpolations between the data points.

Download

In order to understand the origin of this convergence of sea ice volume, we show in Fig. 8 the mean ITDs in March in the Arctic and in September in the Antarctic. March and September correspond to the hemispheric ice extent maxima and are 1 month ahead of the ice volume maxima. The volume gains essentially stop (from S2.09 and S2.07 on in the Northern Hemisphere and Southern Hemisphere, respectively) when there exist categories that can contain the very thick ice produced by ridging (around 10 m and above). When this is the case, the thickest ice occupies only a small fraction of the grid cells, leaving more room for thinner ice that can sustain non-negligible thermodynamic growth rates. Including additional thicker categories allows a more detailed representation of the ice cover but without effects on the growth rates and thereby on the total volumes. The threshold values of 4 m (Arctic) and 2 m (Antarctic) are no coincidence. Maykut and Untersteiner (1971) found that the equilibrium thickness of level ice would reach 3–4 m in standard Arctic conditions. In the Antarctic, a larger climatological snow depth and larger ocean–ice basal heat flux both contribute to a thinner level ice cover on average.

We also note that the Arctic sea ice volume is more sensitive than the Antarctic sea ice volume in the S2 set of experiments (Fig. 7). This is because, in the model, the proportion of Antarctic ice volume produced by deformation (and stored in the thickest categories) is low compared to the Arctic. Creating the categories to host that deformed ice thus has only a moderate effect in the Antarctic. An investigation of the ability of NEMO3.6-LIM3 to realistically simulate the volume of deformed ice is currently underway in another study. We conclude that setting the lower boundary of the thickest category at 4 and 2 m for the Arctic and Antarctic, respectively, is sufficient to allow convergence in simulated ice volumes, provided that the resolution of the ITD in the thin categories remains unchanged. These conclusions hold for present climate conditions and may have to be updated for other applications (e.g., paleoclimate simulations).

3.3 Influence of increasing resolution (S3 experiments)

In the previous set of experiments (S2), we have determined the minimal requirements for the position of the thick ice category in order to achieve convergence of ice volumes, without refining the ITD elsewhere. The final set of sensitivity experiments (S3) is designed to assess the specific role of the resolution of the ITD discretization in the thin range, without changing the upper bound (set to 4 m since we are using a global model). Starting from experiment S2.09 (renamed S3.09 in this third set of experiments), we coarsen and refine the ITD by merging or splitting category boundaries in two. Figure 9 shows the mean seasonal cycles of Arctic and Antarctic sea ice volumes. Increasing the resolution leads to higher sea ice volumes through enhanced bottom growth rates (Fig. A3), as explained in Sect. 3.1. However, the differences become significant in the Arctic only for S3.33, in which the ITD resolution is 4 times finer than in S3.09. Overall, these results suggest that refining the ITD resolution by adding more intermediate categories has a small but not negligible impact on the total simulated sea ice volume.

In line with the results of Sect. 3.2, we conclude that the position of the thickest category (S2 experiments) exerts a first-order control on the total sea ice volume by allowing the existence of thick and deformed ice, while the resolution of the ITD discretization in the thin range (S3 experiments) has a second-order effect on the total volume, by controlling the amount of thermodynamically grown ice.