<figure><img src="image_1.png" /></figure>

HAL
open science

# New representation of sea ice salt dynamics within NEMO-SI3: Evaluation and impacts

Elise Ortega, Martin Vancoppenolle, Clément Rousset, Eva Lemaire

► To cite this version:

Elise Ortega, Martin Vancoppenolle, Clément Rousset, Eva Lemaire. New representation of sea ice salt dynamics within NEMO-SI3: Evaluation and impacts. 2025. hal-05391833

HAL Id: hal-05391833

https://hal.science/hal-05391833v1

Preprint submitted on 1 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

# New Representation of Sea Ice Salt Dynamics within NEMO-SI³: Evaluation and Impacts

Elise Ortega¹ (0009-0000-9628-1326), Martin Vancoppenolle¹ (0000-0002-7573-8582),

Clément Rousset¹ (0000-0002-9722-3969), Eva Lemaire².

¹ Sorbonne Université, CNRS, IRD, MNHN, Laboratoire d'Océanographie et du Climat :

Expérimentations et Approches Numériques, LOCEAN/IPSL, F-75005 Paris, France

² Earth and Life Institute, ELIC, UCLouvain, Belgium

## Key points

13

14

* New options for representing sea ice salt dynamics in NEMO-SI³ were evaluated against a new compilation of ice core salinity observations.

15

16

* Simulated ice salinity variations are consistent with observations, regardless of parameterization details, provided parameter adjustment.

17

18

* Salt dynamics influence sea ice thermodynamics, impacting Arctic sea ice volume by about 10%.
19

**Abstract (250 words)**

Many large-scale sea ice models represent salt dynamics, yet, these remain insufficiently evaluated, with large uncertainty on several parameterization choices. To address this, we conduct global 1° resolution ice-ocean simulations with NEMO-SI³ v5, forced by climatological atmospheric data, and evaluate them against a new compilation of >1,200 ice-core observations. We specifically investigate the effects of (1) the treatment of vertical ice salinity variations (diagnostic vs. prognostic), (2) gravity drainage parameterization type (brine fraction-based vs. Rayleigh number-based), and (3) key parameter values. We show that mean ice salinity can be simulated within 0.1 g/kg (Arctic) and 0.4 g/kg (Antarctic) of the observed mean, provided careful adjustment of parameters. In particular, using a high value for the new ice liquid fraction parameter (75%) proves critical to simulate high salinity in the thin ice range. With appropriately adjusted parameters, the model captures the key large-scale ice salinity features, including the Arctic-Antarctic contrast, seasonal cycle, vertical profiles, and salinity-thickness relationship. A prognostic treatment of vertical salinity increases the variance of ice salinity and ice-ocean salt fluxes, yielding vertically-averaged salinity more in line with observations, however only in the Arctic. The choice of gravity drainage parameterization is influential through interactions with melt ponds, with Rayleigh number-based parameterizations performing best. While salt dynamics strongly affect sea ice mass budget terms, they moderately affect sea ice volume and the ocean due to compensation effects. Uncertainties are larger than in the Antarctic than in Arctic, reflecting less mature process understanding and possibly higher observational biases.
40

## Plain Language Summary

41

42

43

44

45

46

47

48

49

50

51

Sea ice forms from freezing seawater. Unlike freshwater ice, sea ice traps some sea salt when it forms, mostly within tiny pockets of brine. This changes how sea ice grows, melts, and interacts with the ocean. Sea ice models used for climate simulation include these effects, but we do not know well how they capture the way sea ice takes in and releases salt over time. In this study, we look at how on such model represents sea ice salinity. We run simulations to test different ways of simulating salt uptake in and drainage from sea ice and compare the results to a large new dataset of over 1,200 real ice core measurements. We find that models well reproduce observed sea ice salinity patterns, and improve in accuracy when key settings are carefully chosen. We also point out where challenges remain, and give recommendations to modelers on important aspects to consider for improving future simulations.
52

# 1. Introduction

53

54

55

56

57

58

59

60

Sea ice forms from freezing seawater in both polar oceanic regions. Sea ice is a key component of the climate system, influencing air-sea exchanges of heat, mass and momentum (Barry et al., 1993; Chripko et al., 2021; Jenkins & Dai, 2021; Peixoto & Oort, 1992), the ocean circulation (Abernathey et al., 2016; Goosse & Fichefet, 1999; Pellichero et al., 2018; Stössel et al., 2002). Sea ice is also a sensitive indicator of climate change (Stroeve & Notz, 2018). Climate models have represented sea ice processes and their implications for decades (Manabe & Stouffer, 1980; Roach et al., 2020; SIMIP Community, 2020).

61

62

63

64

65

66

67

68

Unlike freshwater ice, sea ice is not only composed of pure ice, but also includes liquid saline brine, hydrated salt minerals and gas bubbles (Hunke et al., 2011; Light et al., 2003; Weeks & Ackley, 1986). This is because sea ice traps some of the salt and seawater constituents upon formation. However, these are hardly incorporated into the ice crystalline lattice. Sea ice salinity, together with temperature, sets sea ice composition, in particular brine fraction (Assur, 1958; Light et al., 2003; Martin Vancoppenolle et al., 2019), which, of all constituents, is best known, in particular regarding effects on sea ice thermal properties (Bitz & Lipscomb, 1999; Notz, 2005; Ono, 1967; Pringle et al., 2007).

69

70

71

72

73

74

75

76

77

78

79

Sea ice salinity exhibits considerable spatiotemporal variability (Cox & Weeks, 1974; Malmgren, 1927). Ice growth mechanisms act as sources of salt, including congelation, frazil and snow-ice growth (Jutras et al., 2016; Notz & Worster, 2009; Smedsrud & Skogseth, 2006) as well as pressure ridge formation (Leppäranta et al., 1995). Conversely, salt sinks are associated with brine drainage (Notz & Worster, 2009), in particular gravity drainage, the convective overturning of brine during the cold season (Notz & Worster, 2009; Wells et al., 2011) and flushing, the percolation of nearly fresh surface meltwater through the brine network during the warm season (Eicken et al., 2002; Untersteiner, 1968; Vancoppenolle et al., 2007). As a result, salinity varies with depth within the ice column and seasonally through time (Eicken, 1998; Malmgren, 1927), with ice thickness (Cox and Weeks, 1974; Kovacs, 1996), and with geographical location (Turner and Hunke, 2015; Vancoppenolle et al., 2009).

80

81

82

83

84

85

The representation of salt dynamics in sea ice and their consequences within the sea ice component of climate models has advanced considerably over the last two decades (Bailey et al., 2020; Boucher et al., 2020; Toyoda et al., 2024) enhancing our understanding of their large-scale effects on ice-ocean exchanges and sea ice thermodynamics. Salt rejection from sea ice and freshwater transport by sea ice drift play a major role in driving water mass transformation at high latitudes, with implications for large-scale oceanic circulation
86 (Abernathey et al., 2016; Goosse & Fichefet, 1999; Pellichero et al., 2018; Stössel et al.,
87 2002). In addition, the temporal variability salt uptake and rejection during sea ice growth
88 has local effects on the surface ocean properties (DuVivier et al., 2021; Uotila et al., 2017;
89 Vancoppenolle et al., 2009). Variations in ice salinity also affect the ice thermal regime,
90 enhancing ice growth and reducing surface melt, with consequences for the seasonal cycle
91 of ice thickness (Bitz & Lipscomb, 1999; Vancoppenolle et al., 2005, 2006). Ultimately, the
92 representation of salt dynamics in sea ice influences mass budget components, the
93 seasonal cycle of sea ice volume, and melt pond coverage (Bailey et al., 2020; DuVivier et
94 al., 2021; Turner & Hunke, 2015; Vancoppenolle et al., 2009).
95
96 Several modeling groups have developed representations of salt dynamics in sea ice,
differing in three respects.

95
96
1. The treatment of vertical ice salinity variations differ. Some models resolve fully
97
98 prognostic salinity profiles, such as CICE/IcePack (Hunke et al., 2024) and MIROC (Toyoda
99 et al., 2024). By contrast, LIM3/SI³ uses a diagnostic profile based on a time-varying mean
100 prognostic salinity (Vancoppenolle et al., 2023). While the fully prognostic approach offers
101 greater vertical detail, whether it provides advantages over a diagnostic profile remains
102 unclear.

103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
2. The parameterization of brine drainage mechanisms, in particular gravity drainage, differ.
With a diagnostic approach to salinity profile, gravity drainage is treated as a restoring term.
With a fully-prognostic salinity profile, more elaborate and physically-based, depth-
dependent formulations are used. In particular, convective parameterizations represent
gravity drainage with a 1D vertical advection equation and a vertical Darcy velocity, for which
the simplest approach is a linear ramp function of brine fraction (Cox & Weeks, 1988, used
in MIROC). More recent parameterizations (Griewank & Notz, 2013; Rees Jones & Worster,
2014; Turner et al., 2013, the latter used in CICE/IcePack) define the Darcy velocity as a
function of a porous-medium Rayleigh number, describing the competition between
destabilizing density forces and viscosity during gravity drainage, and accounting for
permeability (Notz & Worster, 2008; Wells et al., 2011). Alternative, diffusive approaches
also exist (Jeffery et al., 2011; Vancoppenolle et al., 2010), but have been shown to be less
physically realistic for salt and tracer dynamics (Thomas et al., 2020). To date, however, the
performance and large-scale implications of different gravity drainage parameterizations
remain unexplored. (Note there is also some but less diversity in flushing parameterizations,
though they are not considered within the scope of the present study).
3. Calibration of salt uptake and gravity drainage parameters differ, as their choice is not straightforward. The new ice liquid fraction parameter ($\phi_{new}$), which determines the amount of oceanic salt incorporated into newly forming ice, varies across models. It can take a high value, typically 0.75 in CICE/IcePack or much lower values (<0.25) as in LIM3/SI³ and MIROC. More precisely, these latter models adopt a parameterization of the newly formed ice salinity, that assumes instantaneous salt loss upon growth (termed "salt fractionation"), equivalent to using a relatively low $\phi^{new}$ value (Cox & Weeks 1988). However, salt fractionation — and consequently this parameterization of salt uptake — has been invalidated by non-destructive salinity measurements (Notz & Worster, 2009). Recent studies instead favor $\phi^{new} = 0.75$ as physically more appropriate, for various modes of sea ice growth: basal columnar ice, frazil and grease ice, and platelet ice (Smedsrud & Skogseth, 2006; Turner & Hunke, 2015; Wongpan et al., 2021). In addition, gravity drainage parameters are systematically adjusted: parameters regulating drainage intensity, as well as the critical Rayleigh number, vary widely across studies depending on the simulation context (e.g. whether simulations target laboratory or field experiments), simulation length, vertical resolution and time step (Griewank & Notz, 2013; Rees Jones & Worster, 2014; Thomas et al., 2020; Wongpan et al., 2021). In this context, adjusting salt dynamics parameters is needed in large-scale simulations, require careful consideration of both salt uptake and gravity drainage parameters.

To evaluate model parameterizations of salt dynamics, in situ ice salinity observations have been used, though only to a limited extent. One-dimensional modeling studies have employed salinity profile time series, based primarily on ice-core data, but also on impedance-based measurements obtained with a wire harp instrument (Notz et al., 2005), in both laboratory and field experiments (Griewank & Notz, 2013; Rees Jones & Worster, 2014; Turner et al., 2013; Vancoppenolle et al., 2007). Two large-scale modeling studies relied on a literature-based compilation of ice-core salinity observations, consisting of 45 regional monthly mean salinity values, to evaluate the hemispheric and regional seasonal cycle of mean ice salinity (Turner & Hunke, 2015; Vancoppenolle et al., 2009). Finally, several authors have used the salinity-thickness relationships identified in ice-core observational studies (Cox & Weeks, 1974; Kovacs, 1996) to evaluate both one-dimensional and large-scale simulations (Griewank & Notz, 2013; Toyoda et al., 2024; Vancoppenolle et al., 2009). However, available sea ice core salinity observations have not been fully exploited for the evaluation of large-scale salt dynamics, particularly in terms of vertical variations.

In this study, we investigate the representation of salt dynamics in sea ice and their impacts, using newly available vertically resolved parameterizations in the ocean-ice model NEMO.
SI³ v5.0 and a new comprehensive compilation of over 1,200 ice cores (BEPSII, 2025). We analyze large-scale climatological 1° ice-ocean simulations, with differing representations of salt processes. In particular, we examine (1) the vertical treatment of ice salinity (diagnostic vs. prognostic), (2) the gravity drainage parameterization type (brine fraction-based vs. Rayleigh number-based), and (3) key parameter values. We evaluate the large-scale ice salinity features, spatially (both horizontally and vertically), seasonally, and through the salinity-thickness relationship. We also assess the impacts on sea ice mass budget, volume and ice-ocean salt exchanges. Our ultimate objective is to identify the optimal configuration for NEMO-SI3, while providing insights of broader relevance to the sea ice modeling community.

We detail the methods used in Section 2, evaluate simulations against observations in Section 3.1, and assess their impact on sea ice thermodynamics in Section 3.2, and on ice-ocean exchanges in Section 3.3. The results are discussed in Section 4.2.

## 2. Materials and methods

### 2.1 The ocean sea ice model NEMO-SI³

We run simulations with the ocean modelling system NEMO (Nucleus for European Modelling of the Ocean) version 5.0 (Madec et al., 2022). The ocean component is a three-dimensional, free-surface, hydrostatic, Boussinesq, primitive equation model, solved with finite-difference numerics. A nonlinear free surface approach (z* coordinate) allows cell thicknesses throughout the water column to vary with time, which favors an accurate representation of the surface freshwater flux (Adcroft et al., 2004).

The sea ice component is SI³ (Sea Ice modeling Integrated Initiative, Vancoppenolle et al., 2023), a dynamic-thermodynamic continuum model, with multi-layer energy and salt conserving snow and ice thermodynamics, the adaptive elastic-viscous-plastic (aEVP) rheology for ice dynamics (Kimmritz et al., 2016), a second order moment conserving scheme for horizontal advection (Prather, 1986), and a multi-category framework to represent subgrid-scale variations in ice thickness (Bitz et al., 2001; Lipscomb, 2001; Thorndike et al., 1975).

### 2.2. Sea ice salt dynamics and their impact

Sea ice is assumed to consist of liquid brine and solid ice, with all the salt contained in the liquid brine. Temperature (T) and salinity (S) characterize the thermodynamic state of sea
ice. Under the thermal equilibrium assumption, brine fraction is given by $\phi(z) =$ $S/S_{br}(T)$, with $S_{br}$ the brine salinity, which follows a third-order liquidus relationship (Vancoppenolle et al., 2019). Like temperature, sea ice salinity evolves in time and space, and is defined for all categories and vertical layers. Processes such as horizontal transport, ridging, and rafting — which shape variations in physical and thickness space — are considered. For these processes, salinity is converted into its extensive, conservative counterpart, i.e., salt content. Vertical and temporal variations in ice salinity mostly arise from salt uptake and drainage, which are treated as sources and sinks of salt.

## 2.2.1 Vertical salinity variations

At this stage, it is worth mentioning the two options available in SI³ for representing vertical salinity variations. In the fully prognostic salinity profile approach — introduced in NEMO v5.0 — depth-dependent formulations of brine drainage are used to calculate temporal variations in ice salinity at each depth. In contrast, the diagnostic profile approach (used as the reference in earlier NEMO versions), computes variations in vertically averaged salinity ($\hat{S}$) using empirical formulations of brine drainage, while the vertical salinity profile is diagnosed from $\hat{S}$. Specifically, salinity is assumed vertically uniform when $\hat{S}$ exceeds 4.5 g/kg, linear from 0 g/kg at the surface to $2\hat{S}$ at the bottom when $\hat{S}$ is below 3.5 g/kg, and defined as a combination of both cases in the intermediate range.

## 2.2.2. Salt uptake

The incorporation of salt into sea ice is tied to sea ice growth, with contributions from all sources of sea ice mass: basal growth, new ice formation in open water, snow ice, and the entrapment of seawater into newly formed ridges. For each growth process x, the salinity of the new ice is formulated as

$$S_{new}^{x} = \phi_{new}^{x} \cdot S_{w},$$

where $S_w$ is seawater salinity of and $\phi_{new}^x$ is the liquid fraction of newly formed ice. The salt content of the newly formed ice is conservatively combined with that of pre-existing ice.

Under the diagnostic salinity profile approach, this modifies the vertically averaged salinity.

The value of $\phi_{new}^x$ depends on the specific ice growth process. For snow-ice formation, $\phi_{new}^x$ corresponds to the air volume initially present in the snow, diagnosed from the density difference between ice and snow (Jutras et al., 2016). For new ridges, $\phi_{new}^x = 0.3$ corresponds to the prescribed new ridge porosity parameter (Leppäranta et al., 1995).
For new ice formation in open water and at the ice base, several options exist. Before the version 5 of NEMO, $\phi_{new}^x$ was calculated as a function of basal growth rate following (Cox & Weeks, 1988); and for open water growth, as a function of new ice thickness (Kovacs, 1996). These formulations yield relatively low values of $\phi_{new}^{bg}$ around 0.1 for basal growth and $\phi_{new}^{ow} = 0.4$ for new ice formation in open water. However, several studies suggest that a value of 0.75 better represents the physics of new ice growth across frazil/grease ice, columnar and platelet growth modes (Smedsrud & Skogseth, 2006; Turner et al., 2013; Wongpan et al., 2021). Consequently, we now set $\phi_{new}^x = 0.75$ for both basal and open water growth modes. The implications of this choice are evaluated through targeted sensitivity experiments (see Section 2.4).

### 2.2.3. Brine drainage

The formulation of gravity drainage depends on the selected approach for representing vertical salinity variations. Table 1 details the gravity drainage parameterizations available in NEMOv5.0 and evaluated in this work.

In the diagnostic profile approach, both gravity drainage (and flushing) are parameterized using exponential restoring terms. The asymptotic salinities are set to 5 g/kg for gravity drainage and 2 g/kg for flushing, whereas the restoring time scales are set to 20 days for gravity drainage, and 10 days for flushing, following values established in Vancoppenolle et al. (2009)

In the prognostic salinity profile approach, the time evolution of the vertical sea ice salinity profile is governed by a 1D advection equation:

$$ \frac{\partial s}{\partial t} = -(w_{gd} + w_{flu}) \frac{\partial s_{br}}{\partial z} $$

where $w_{gd}$ and $w_{flu}$ are the vertical Darcy velocities (positive upwards) associated with gravity drainage and flushing, respectively. Three gravity drainage parameterizations were implemented in SI³ for NEMOv5.0. Following Thomas et al. (2020), only convective parameterizations of gravity drainage were considered. The first parameterization (Cox & Weeks, 1988; as reformulated by Thomas et al., 2020), the gravity drainage Darcy velocity is defined as a linear ramp function of liquid fraction. The two other retained convective gravity drainage parameterizations (Griewank & Notz, 2013; Rees Jones & Worster, 2014) are physically more elaborated, as they specify the convective velocity as a function of a local porous-medium Rayleigh number (Notz & Worster, 2009). The Rayleigh number, Ra, characterizes the competition between destabilizing density forces and viscosity during gravity drainage, accounting for permeability, and is defined as:
$$249 \qquad Ra(z) = \frac{cg\beta [S_{br}(z) - S_w]}{kv} \bar{\bar{I}} (h-z).$$

c, k, v are the heat capacity, thermal conductivity, and kinematic viscosity of brine, respectively; g is the acceleration due to gravity, β is the haline density coefficient, and $\bar{\bar{I}}$ is the ice effective permeability. The two Rayleigh-based parameterizations differ in the functional dependence of the Darcy velocity on $Ra$ in the convecting layer, as highlighted by Worster and Rees Jones (2015). In the Griewank & Notz (2013) parameterization, the Darcy velocity is proportional to the integral of the local Rayleigh number across the convective region, whereas in the Rees Jones & Worster (2014) parameterization, it is proportional to the maximum value.

Flushing is represented following the assumption that ice is permeable to percolation of surface meltwater, if any available, provided the liquid fraction is above $\phi_c = 5\%$ throughout the entire ice thickness (Vancoppenolle et al., 2007). A fraction $\alpha_{flu} = 30\%$ of the surface meltwater flux $F_{fw}^{su}$ (calibrated against observations of early summer salinity profiles) percolates instantaneously through the ice, displacing the brine network downwards. The corresponding Darcy velocity can be formulated as $w_{flu} = -\alpha_{flu} F_{fw}^{su}/\rho_{br}$, or zero if $\phi < \phi_c$ in any of the ice layers.

To solve the 1D vertical salt advection equation, a time-splitting approach with sub-iterations is used. The number of sub-iterations is determined as a function of the maximum Courant–Friedrichs–Lewy (CFL) condition across all vertical layers, defined as $CFL_k = |w_k| \times \Delta t_{ice}/\Delta z_k$, where $w_k$ is the brine Darcy velocity, $\Delta z_k$ the thickness of the kth layer, and $\Delta t_{ice}$ is the sea ice model time step. To prevent endless loops, the number of sub-iterations is capped at 50 for gravity drainage and 10 for flushing. The required number of sub-iterations increases with vertical resolution. Numerical instabilities may therefore arise when a very large number of vertical layers are used. In the simulations presented here, which use 10 vertical sea ice layers, we never needed more than 10 iterations for gravity drainage and 2 for flushing.

### 2.2.4. Impacts of salt dynamics

In SI³, sea ice salinity influences temperature, growth and melt, and melt ponds. Ice salinity influences the vertical temperature profile within sea ice because salinity affects the sea ice thermal properties. This can be seen in the heat diffusion equation, from which the temporal evolution of the vertical temperature profile is calculated:

$$280 \qquad \rho c_p(S, T) \frac{\partial T}{\partial t} = \frac{\partial}{\partial z} \left[ k(S, T) \frac{\partial T}{\partial z} \right] + Q_{sw}^{abs},$$
where, $Q_{sw}^{abs}$ is the absorbed solar radiation per unit volume, $\rho$ is density and both the specific heat $c_p(S, T)$ and thermal conductivity $k(S, T)$ are formulated as functions of $S$ and $T$, in order to account for the effects of brine (Ono, 1967; Pringle et al., 2007). $c_p(S, T)$ largely increases near the freezing temperature, whereas salinity implies a mild decrease in thermal conductivity, as compared with pure ice values.

Ice salinity also affects sea ice growth and melt rates. This occurs because, in SI³, the sea ice enthalpy $q$ — which determines the energetic cost of melting and the heat release during freezing — is a function of $S$ and $T$, accounting for brine effects on the latent and sensible heat required to form or melt sea ice. Indeed, the mass per unit area added to or removed from sea ice $\Delta M$ (kg/m²) is given by Schmidt et al., 2004:

$$ \Delta M = \delta Q \Delta t / \Delta q. $$

Here $\delta Q$ is the net heat flux forcing for each freezing / melting process applied during time step $\Delta t$ and $\Delta q$ is the enthalpy difference involved in the phase transition — either freezing or melting. Ice salinity affects the enthalpy of forming or melting ice, increasing or reducing the sea ice mass change, for a given forcing.

Finally, we use a representation of melt ponds, specific to SI³, combining elements from previously presented melt pond parameterizations (Flocco et al., 2010; Holland et al., 2012; Hunke et al., 2013) including pond growth and drainage, but neglecting lateral drainage and topography effects (Sterlin et al., 2021). In each category, a melt pond reservoir can store some melt water, be covered by a frozen lid, and is advected and redistributed as other conservative sea ice tracers. The volume of melt water in ponds per unit grid cell area ($V_p$), evolves according to runoff, overflow, drainage and freezing. Runoff is a source term including the contributions of surface snow, ice and lid melting. Overflow is a sink term formulated such as pond depth and fraction cannot exceed threshold values. Drainage is parameterized following Flocco et al., 2010, and depends on sea ice permeability. Following Holland et al. (2012), pond area fraction and depth are derived from pond volume, assuming that pond depth linearly scales with pond area. Melt ponds largely reduce surface albedo as compared with a bare-ice case, in proportion with pond fraction (Sterlin et al., 2021). Overall, the representation of ice salinity can affect melt ponds through its impact on surface melt rate and permeability (affecting drainage).

## 2.3. Ice salinity observations

As an observational reference, we use a new ice salinity and temperature ice core database (BEPSII, 2025). With 1,236 ice cores from both hemispheres, this represents the largest ice
salinity database available to date. We exclude ice cores with more than 60% missing
depths or with total section length shorter than 60% of ice thickness, leaving 1,113 ice cores
— 391 from the Northern hemisphere and 768 from the Southern hemisphere (Figure 1). Ice
cores with salinity exceeding 2 g/kg in the uppermost 20 cm typically exhibit C-shape salinity
profile characteristic of the action of gravity drainage, hence we classify them as high-
surface-salinity (HSS, N = 954). By contrast, low-surface salinity cores (LSS, N = 159)
indicate the action of flushing, which is particularly prevalent during the Arctic summer. This
classification was adjusted and verified visually, on a profile-by-profile basis. The salinity-
thickness relationship is evident in the 948 HSS cores (as defined above), yielding <S> =
5.14 + 0.52 / h (R² = 0.31), which closely matches the findings of Kovacs (1996), who
reported S = 4.90 + 0.91 / h (R² = 0.73), based on several hundreds of first-year ice cores.

The ice core database contains several biases. The well-known measurement bias due to
brine loss upon core extraction results in an underestimation by up to 5 g/kg in the bottom 10
cm of the ice cores (Eicken et al., 2010; Notz et al., 2005). There is also a hemispheric
sampling bias, Antarctic ice cores accounting for two thirds of the total database.
Additionally, the distribution of ice thickness is not uniform, with 85% of the cores having
lengths ranging from 0.2 m to 2 m, leaving very few cores with either short (<0.2 m) or
relatively long (> 2m) core lengths. Finally, a seasonal sampling bias also exists, with, for
instance, only 6 samples collected in October in the Arctic compared to 236 collected in
November in the Antarctic.

## 2.4. Experimental setup and tuning

A series of NEMO-SI³ simulations with varying salt dynamics parameterizations were
conducted, five of which are retained in the final analysis (Table 2). These simulations use a
global eORCA1 tripolar isotropic grid, with nominal 1° horizontal resolution (refined to 1/3° at
the equator) and 75 vertical levels increasing from 1 m depth near the surface to 200 m at
6000 m depth. They are forced by the CORE normal year atmospheric dataset (Large &
Yeager, 2004), and initialized from a 1000-year spin-up simulation from de Lavergne et al.
(2025). Sea surface salinity restoring is activated everywhere except beneath sea ice. We
use 10 ice categories, 5 vertical snow layers and 10 ice layers, and set a minimum ice
thickness to 5 cm. Simulations are 10-year long, which is sufficient for sea-ice to reach
equilibrium. In most cases, year 10 is analyzed. All simulations exhibit the typical low
Antarctic sea ice extent bias characteristic of forced-atmosphere simulations, with minimum
sea ice extent values of around 1 million sq km instead of the 3-4 million sq km expected
from observations.
348 The simulations focus on accurately representing salt dynamic, with an emphasis on the new
349 ice liquid fraction parameter $\phi_{new}$, the handling of vertical salinity variations, comparing
350 prognostic and diagnostic profiles, and the selection of gravity drainage parameterization.
351 For each simulation, the gravity drainage parameters were fine-tuned based on the observed
352 salinity-thickness relation-ship (Figure 2), aiming to closely match the asymptotic salinity
353 value of 5.1 ± 0.1 g/kg for cold first-year thick ice.

354 In the 5 simulations studied, RESTO_REF and RESTO_MOD use diagnostic vertical salinity
355 profiles while CW, GN and RJW use prognostic profiles. RESTO_REF uses a low value of
356 $\phi_{new}$= 0.25, consistent with the reference setup preceding version 5 of NEMO, whereas
357 RESTO_MOD increases $\phi_{new}$ to 0.75, considered more physically realistic. In
358 RESTO_MOD, the gravity drainage timescale was reduced by a factor of 10 to achieve the
359 same thick ice salinity as NEMO_REF. The CW, GN and RJW simulations incorporate
360 vertically-resolved salt dynamics, each paired with a distinct convective gravity drainage
361 parameterization (Table 1). For these simulations the intensity of brine convection ($\alpha$) and
362 the critical Rayleigh number ($Ra_c$) were tuned. The resulting values of $\alpha$ and $Ra_c$ differ from
363 those in Thomas et al., 2020, whose tuning was based on short-term (< 10-day) one-
364 dimensional simulations designed to be match artificial young ice time series observations.
365 In contrast, our tuning yields weaker and more intermittent gravity drainage, aligning better
366 with yearlong simulations of natural Antarctic sea ice (Wongpan et al., 2021).

367 **3. Results**

368 **3.1. Salt Dynamics**

369 **3.1.1 Ice Salinity-Thickness Relationship**

370 We begin by analyzing the salinity-thickness relationship in both hemispheres, comparing
371 simulations to observations, as this relationship guided the calibration of ice salinity.
372 Simulated salinities, vertically averaged and resolved by thickness category, are presented
373 as 2D histograms (Figure 2), alongside ice core salinities plotted against core length. To
374 focus on cold growing ice, we excluded the categories with low simulated surface salinity (<2
375 g/kg) and the LSS observations from the analysis.

376 A notable feature is that in RESTO_REF, the vertically-averaged salinity remains below 9
377 g/kg. In contrast, the other simulations, consistent with observations, show salinity values
378 significantly exceeding 9 g/kg for ice thinner than 1 m, with peaks reaching ~25 g/kg for ice
379 around 10 cm thick. This discrepancy is attributed to the relatively low value of $\phi_{new}$ in
380 RESTO_REF (0.25), which limits salt uptake compared to the other simulations (all using
381 $\phi_{new}$ = 0.75).

382 A second important aspect is that the inverse relationship between salinity ($S$) and thickness
383 ($h$), as identified in the observations (see Section 2.3), is also reproduced in the simulations.
384 For each simulation, an $S-h$ curve is fitted using a function of the form $A+B/h$ (Figure 2). The
385 variance explained by this fit exceeds 75% in most simulations ($R^2$ = 0.75-0.85), though it is
386 lower in RESTO_REF ($R^2$=0.57). In all cases, the variance remains higher than in the
387 observations ($R^2$ = 0.31), which may reflect observational limitations such as insufficient
388 sampling, random measurement errors, or strong natural variability not fully captured by the
389 simulations.

390 The asymptotic thick ice salinity value $A$ lies within 0.1 g/kg of the observation-based
391 reference (5.13 g/kg) for all simulations except for RESTO_MOD where it is slightly
392 underestimated (4.18 g/kg). The narrow range in the values of $A$ reflects the effort on tuning
393 the salinity of cold, thick ice. In contrast, the slope coefficient $B$ shows greater variability:
394 observations yield a lower value (0.52 g/kg.m), as does RESTO_REF (0.33 g/kg.m),
395 compared to the other simulations (0.90-97 g/kg.m). These differences in $B$ stem from a
396 spread in ice salinity for ice thinner than 1 m, where simulated salinity depends on the
397 balance between salt uptake and gravity drainage, both of which differ among simulations.
398 Additionally, observed salinity for thin ice is subject to more uncertainty than for thick ice.
399 Thin ice is warmer, saltier, and more permeable, making it more prone to salt loss during
400 core extraction and complicating comparisons with simulations.

401 Further differences among simulations are evident in the density plots of Figure 2.
402 Simulations using a prognostic salinity profile and convective gravity drainage (CW, GN,
403 RJW) show higher salinity variance around the $S-h$ curve, aligning more closely with ice core
404 observations than simulations with a diagnostic profile and restoring gravity drainage
405 (RESTO_REF, RESTO_MOD). This distinction is especially pronounced for ice thicker than
406 1 m.

### 3.1.2. Hemispheric and seasonal salinity variations

408 Next, we investigate the mean values (Table 3) and seasonal cycles (Figure 3) of the
409 vertically averaged salinity, for each hemisphere. Overall, all the simulations capture
410 hemispheric mean salinities reasonably well. In the Northern hemisphere, RESTO_MOD,
411 GN, and RJW align closely with the observed mean salinity (4.4 g/kg), deviating less than
412 0.1 g/kg. In contrast, RESTO_REF underestimates it (3.8 g/kg) while CW overestimates it
413 (5.1 g/kg). In the Southern hemisphere, CW (8.0 g/kg), GN (7.8 g/kg), and RJW (7.6 g/kg)
414 significantly exceed the observed value (6.6 g/kg). Conversely, RESTO_REF (5.3 g/kg) and,
415 to a lesser degree, RESTO_MOD (6.2 g/kg) underestimate it.

416 All simulations, consistent with observations, predict higher ice salinity in the Antarctic
417 compared to the Arctic. The Antarctic-Arctic mean ice salinity difference ranges from 1.5 to
418 3.3 g/kg across simulations (compared to 2.2 g/kg in observations). This difference arises
419 from several mechanisms captured in all simulations: flooding and snow-ice formation are
420 frequent in the Antarctic but rare in the Arctic, whereas low-salinity multi-year ice is
421 widespread in the Arctic and largely absent in the Antarctic.

422 In both hemispheres, the simulated seasonal cycle of sea ice salinity is also overall
423 consistent with observations (Figure 3). All simulations reproduce a salinity maximum in fall,
424 associated with the formation of large amounts of new ice that traps salt, and a summer
425 minimum, resulting from extensive salt loss through brine drainage, particularly via flushing
426 in the Arctic. In addition, the simulated salinities fall within the observed range for most
427 months, with a few exceptions. In particular, simulated ice salinity is higher than observed
428 during the Arctic summer, and lower during the early growth phase (October in the Arctic
429 and March in the Antarctic). However, monthly differences should be interpreted with
430 caution, as the largest discrepancies coincide with periods of lowest observational coverage.
431 A notable distinction among simulations is that, during the Antarctic cold season,
432 RESTO_MOD produces salinities that are lower and more consistent with observations than
433 CW, GN and RJW. Finally, the tendency for RESTO_REF and CW to simulate the lowest
434 and highest salinities, respectively, persist throughout most months.

435 Table 3 highlights differences in salinity variance across simulations. In the Arctic,
436 simulations using a prognostic salinity profile (CN, GN, RJW) exhibit a mean salinity
437 variance (1.7 g/kg) close to observations (1.8 g/kg), whereas those using a diagnostic
438 salinity profile show lower variance, particularly RESTO_REF (0.9 g/kg). In the Antarctic, all
439 simulations reproduce the observed salinity variance (2.7 g/kg) within 0.2 g/kg, except
440 RESTO_REF, which shows a lower value (1.0 g/kg).

### 3.1.3. Vertical salinity profiles

442 We now examine vertical salinity variations. The high surface salinity (HSS) and low surface
443 salinity (LSS) profile types identified in observations are also reproduced in the simulations.
444 The four panels of Figure 4 show the average normalized profiles for both surface salinity
445 classes in each hemisphere. The simulated frequencies of each class are commensurate
446 among simulations and observations (Table 3): 35-54% of profiles exhibit low-salinity in the
447 Northern hemisphere (compared to 35% in observations) while 2-7% do so in the Southern
448 hemisphere (compared to 7% in observations). We do not expect perfect agreement with
449 observations for this metric, as our classification criterion is sensitive to profile shape, which
450 varies substantially between simulations using diagnostic and prognostic approaches.

451 Based on their ability to reproduce vertical profile shapes, the simulations can be divided in
452 two groups. As expected, those using a prognostic salinity profile (CW, RJW, GN) resolve
453 ice salinity variations in greater vertical detail compared to those using a diagnostic salinity
454 profile (RESTO_REF, RESTO_MOD), which represent salinity in a more idealized manner.
455 In CW, RJW and GN, the HSS class shows a C-shaped profile, characterized by high salinity
456 near the interfaces, and lower salinity in the ice interior. The LSS class closely resembles a
457 typical summer profile, featuring nearly fresh ice near the surface (0.5–1 g/kg) and more
458 saline ice at depth (4–6 g/kg). In contrast, the RESTO simulations also reproduce a winter-
459 summer contrast in profile shape, but with less detail: HSS profiles are vertically uniform,
460 while LSS profiles are linear, as prescribed a priori by the diagnostic method.

461 A striking difference evident in Figure 4 is that prognostic parameterizations exhibit
462 significantly higher bottom salinity than observed in simulations. In the Arctic, it is well
463 established that the low salinity bias caused by brine loss during core extraction is largest
464 near the ice base (see Section 2.3; Notz et al., 2005, and Section 2.3). Therefore, Arctic
465 simulated salinity profiles may be more realistic in CW, GN and RJW than Figure 4
466 suggests. In the Antarctic, the model-observation discrepancy at the ice base is even larger
467 than in the Arctic. However, our understanding of ice desalination physics and salinity
468 measurements biases remains less developed. Finally, we note that diagnostic
469 parameterizations cannot exhibit this bottom-intensified salinity difference, as they prescribe
470 vertical salinity profiles a priori.

471 Observations suggest that the ice salinity profile shape slightly differs between the Arctic and
472 the Antarctic, and only simulations using a prognostic salinity profile reproduce this. More
473 specifically, CW, GN and RJW simulate HSS profiles with higher near-surface salinity values
474 in the Antarctic (9–12 g/kg) than in the Arctic (~7 g/kg), a tendency most pronounced in
475 RJW. Also, these simulations reproduce the greater non-linearity of LSS profiles in the
476 Antarctic, where they show a near-surface inflection, in contrast to the more linear profiles
477 typically observed in the Arctic.

### 3.1.4. Geographical distribution

479 The geographical distribution of vertically averaged salinity (Figure 5) is broadly consistent
480 across simulations and aligns with previous modeling studies (Turner & Hunke, 2015;
481 Vancoppenolle et al., 2009). GN and RJW behave very similarly, and CW shows only
482 localized differences from these two. In contrast, RESTO_REF and RESTO_MOD differ
483 more substantially. However, a rigorous observation-based assessment of salinity maps is
484 hindered by limited spatial sampling in the observational database.

485 During the Northern hemisphere winter, all simulations depict a contrast between relatively
486 saline seasonal ice (S > 5 g/kg) along the Arctic periphery and fresher perennial ice (S <
487 5 g/kg) within the Arctic Basin. However, differences emerge in the magnitude of salinity
488 values: diagnostic parameterizations show weaker spatial contrasts, with RESTO_MOD
489 showing the most homogeneous ice due to its lower gravity drainage time constant and
490 consequently stronger drainage. By contrast, prognostic parameterizations produce greater
491 spatial variability. Regionally, GN and RJW simulate less saline ice than CW in Hudson Bay,
492 the Laptev Sea, and the Central Arctic. In summer, Arctic ice salinity is more uniform, though
493 a small patch of higher salinity off West Siberia appears with prognostic parameterizations,
494 especially in CW. Simulated salinities in RESTO_REF and RESTO_MOD remain relatively
495 low (≈2 g/kg) and homogeneous across the basin.

496 In the Antarctic, prognostic parameterizations exhibit higher and more variable ice salinity
497 than the diagnostic ones. Since most ice melts in summer, we only show maps at maximum
498 extent. The largest inter-simulation differences occur in the outer belt of thin (<1 m), young
499 ice, where CW, GN, and RJW reach 6–10 g/kg, while RESTO_REF ranges from 5–8 g/kg,
500 and RESTO_MOD from 5–6 g/kg. Among these, GN and RJW produce the highest salinities
501 within the outer thin ice belt, and show a stronger contrast with thick ice in the western
502 Weddell and Ross Seas compared to CW. Due to its low gravity drainage time constant,
503 RESTO_MOD features a narrow circumpolar belt of saline young ice, while RESTO_REF,
504 with slower drainage, displays a wider belt.
## 3.2. Impacts of salt dynamics

To assess the large-scale effects of sea ice salinity, the hemisphere-integrated ice volume is calculated for each simulation. Results indicate that simulations with more saline ice feature an overall larger ice volume, in both hemispheres (Figure 6). In the Northern hemisphere (NH), the annual mean total ice volume varies by ~2,100 km³ across simulations.

RESTO_REF yields the lowest volume, while CW produces the highest, and RESTO_MOD, GN, and RJW fall in between. These differences remain relatively stable throughout the year. In the SH, volume differences are smaller (<1,000 km³) at the time of maximum volume (October). RJW, GN and CW have higher volumes than RESTO_MOD and RESTO_REF. The greater volume differences across simulations in the NH reflect the presence of multi-year ice in the Arctic, where differences can accumulate over several years.

Inter-simulation volume differences stem from a complex interplay between sea ice mass and salt budget terms, accumulated over the 10 years of simulation. During this time, the sea ice mean state adjusts so that total annual melt balances growth. Classically, under constant salinity conditions, sea ice adjusts its thickness distribution to achieve mass budget equilibrium (Bitz & Roe, 2004). However, with interactive salinity, both thickness and salinity are involved in the adjustment process, adding an additional layer of dynamical feedback.

To better understand ice volume differences between simulations, the annually integrated mass budget is diagnosed and analyzed. Each term of the mass budget (bottom, snow-ice growth, and open water ice growth, bottom melt, surface melt and lateral melt) is integrated annually, for each hemisphere (Figure 7). By year 10, the net annual mass budget imbalance, averaged across all simulations, reaches 0.04 (0.28) x 10³ Gt ice / year in the NH (SH), corresponding to 0.1% (2.7%) of the total hemispheric ice mass (Figure 7 a-c).

Similarly, the salt budget imbalance amounts 0.7% of the total hemispheric salt mass (Figure 7 b-d). Therefore, regarding mass and salt balance, our simulations can be considered close to equilibrium.

### 3.2.1 Impacts on sea ice mass budget

Differences in the mass budget terms across simulations can be both positive and negative, and are primarily driven by the thermodynamic effects of salt dynamics. When integrated over the NH (SH), the total growth and melt vary across simulations by 1.0 (0.9) x 10³ Gt of ice, equivalent to 5% of the total. However, individual mass budget terms, especially open water ice growth, bottom growth and bottom melt, exhibit greater variability, differing by up to 20% of total annual growth. Due to compensating positive and negative contributions, these
changes tend to offset each other in the net budget. Regarding the drivers of the changes in
sea ice mass budget, the ocean's role appears negligible, as differences in sensible ocean-
to-ice heat flux across simulations remain mostly within 0.1 W/m². Instead, changes in
individual mass budget terms seem essentially driven by direct effects of salinity on thermal
properties (q, c_p, k). In some cases, indirect feedbacks, such as the growth-thickness
relationship and the influence of melt ponds on surface albedo, also contribute to these
changes.

The largest change in mass budget is due to $\phi_{new}$. Annual freezing (i.e., the sum of bottom, open water, and snow ice growth) increases by 1.6 (2.3) x 10³ Gt in the NH (SH) when $\phi_{new}$ increases from 0.25 to 0.75. This effect is isolated by comparing RESTO_MOD to RESTO_REF. The difference reflects a more than twofold increase in open water growth — by 3.9 (4.6) x 10³ Gt of ice — and a concurrent reduction in bottom growth — by 2.3 (2.4) x 10³ Gt of ice. In RESTO_MOD the larger $\phi_{new}$ reduces the enthalpy of new ice, lowering the thermal cost of new ice formation by a factor of ~3, which enhances open water growth. The same enthalpy-mediated effect applies to bottom growth. However, the simultaneous reduction in thermal conductivity dominates, resulting in a net decrease in bottom growth.
Finally, basal melting mitigates the effect of anomalous growth on ice volume, as ice melts more in RESTO_MOD than in RESTO_REF, by 1.1 (1.6) x 10³ Gt of ice in the NH (SH), which we attribute to the lower enthalpy of the more saline ice in RESTO_MOD.

Using a prognostic instead of a diagnostic profile intensifies basal ice growth and melt, though this effect — highlighted when comparing the prognostic parameterizations to RESTO_MOD — is smaller than the impact of $\phi_{new}$. The causes are not obvious to identify.
First, open water growth weakly differs (within 0.3 x 10³ Gt or ice) across these simulations, as it is primarily controlled by $\phi_{new}$, which remains identical in all cases. By contrast, bottom growth increases more substantially, by 0.8 (1.30) x 10³ Gt of ice, when a prognostic profile is used. In the Arctic, this increase relates to a lower salinity in the Arctic Basin in the prognostic parameterizations (by 0.5-1 g/kg), which increases the ice thermal conductivity compared to RESTO_MOD. In the Antarctic, the cause of increased bottom growth is less straightforward, though the most plausible explanation is a more stable ocean when using prognostic salinity, potentially due to reduced salt rejection. Additionally, bottom melt is higher with a prognostic profile, by 0.7 (1.20) x 10³ Gt of ice, consistent with the presence of more saline ice at the base when a prognostic profile is resolved.

The effects of the gravity drainage parameterization mass budget are generally small. For instance, GN and RJW remain extremely similar in terms of volume and mass budget terms — with differences never exceeding 0.1 x 10³ Gt of ice. A notable exception is CW, which
simulates larger Arctic sea ice volume than RJW and GN, by 2,100 km³ (Figure 6), the largest inter-simulation volume difference, despite being comparable to GN and RJW in terms of mass budget terms. This suggests a unique equilibration process in CW, emphasizing the importance of the transition to equilibrium. To investigate this further, we analyze the first year of the simulation (Figure 8) which reveals that CW develops large Arctic sea ice volume differences compared to other simulations as early as the first melt season.

It appears that the choice of the gravity drainage parameterization can largely impact the Arctic sea ice volume through melt pond feedbacks. In CW, the low melt pond coverage plays a central role in a series of feedbacks that shape how sea ice volume equilibrates. The origin of these differences lies in CW producing more saline near-surface ice (by 1–2 g/kg, Figure 8b) at the end of the growth season. The higher salinity increases the thermal inertia of the surface, making it more resistant to warming. Consequently, surface ice in CW remains colder by 0.4–0.8°C during the early melt season (Figure 8c). This colder ice surface enhances heat conduction and sea ice growth at the end of the growth season, leading to thicker ice (Figure 8d). During the melt season, thickness differences amplify because the colder ice in CW delays surface melt (Figure 8e), resulting in much lower melt pond fraction compared to other simulations (Figure 8f). This, in turn, reduces surface melt. Note that as surface melting and melt pond formation remain extremely rare in the Southern hemisphere, this feedback process does not occur here, and CW remains closer to other simulations in that region.
### 3.2.2 Impacts on the ice-ocean exchanges

While previous works suggest a weak impact of the representation of salt dynamics on the ocean (see, e.g. Uotila et al., 2017), we further investigate the impact of the representation of salt dynamics on the ocean. To this end the sea ice salt budget and its contributors was diagnosed and analyzed. The net ice-ocean salt flux was decomposed into ice growth, melt and brine drainage contributors, then annually integrated, for both hemispheres and each simulation (Figure 7 b-d).

We find that while the gross contributions to the ocean-ice salt flux largely vary across simulations, the net flux is much more consistent. Indeed, the gross salt flux budget terms are much lower in RESTO_REF than in the other simulations. A key difference is the markedly reduced salt uptake during ice growth, which is on average 36 % (39 %) lower in the NH (SH). This reduction primarily results from the lower value of $\phi_{new}$, and to a lesser extent from reduced overall ice growth in RESTO_REF. Consequently, less salt is available for release through brine drainage and melting. By comparison, net salt gains and losses remain relatively consistent among simulations. In RESTO_REF, the net annual salt gain reaches 90 % (79 %) of the average value in the other simulations, as drainage and melt processes tend to offset uptake contributions on timescales of a few days.

The net ice-ocean salt flux — including all growth and melt processes as well as brine drainage — roughly scales with the ocean-ice volume flux across all simulations (Figure 9, top five panels). However, the representation of ice salinity clearly affects this relationship: a prognostic salinity profile shows larger spread in the salt versus volume flux compared to a diagnostic approach. Despite these differences, the overall impact on the ocean remains limited. This is because changes in ocean salinity are strongly influenced by the mass fluxes associated with ice-ocean exchanges, which occur alongside salt fluxes. To illustrate this, we calculate the salinity tendency of an upper ocean layer of thickness $h_w = 20$ m and of salinity $S_w = 34$ g/kg, as:

$$ \frac{dS_w}{dt} = 1/h_w \cdot [F^m S_w - F^s / \rho_w \times 1000] , $$

where $F^m$ and $F^s$ are ice-ocean mass and salt fluxes, respectively, and $\rho_w$ is seawater density. The bottom panels of Figure 9 show that the net desalination rate of the upper ocean scales quasi-linearly with the mass flux, with little variation in this relationship across simulations. This highlights the relatively small role of salt fluxes compared to mass fluxes in setting the surface ocean salinity tendency.
626

# 4. Discussion and conclusions

627

628

629

630

631

632

633

634

We evaluate newly implemented, vertically-resolved salt dynamics in sea ice in NEMO-SI³, a state-of-the-art ocean-sea ice model. A key novel aspect is the comprehensive evaluation of two Rayleigh-number based gravity drainage parameterizations, against observational data and an established reference approach. These parameterizations are tested here for the first time in a large-scale setup, using the same model, configuration and tuning targets.

Moreover, our observational reference is a new ice salinity database based on >1,200 cores, largely surpassing previous efforts.

Our analysis highlights key elements in the representation of salt dynamics in sea ice.

635

636

637

638

639

1. Calibrating gravity drainage parameters is essential to adjust simulated mean salinity, and the salinity-thickness (*S*-h) relationship is a valuable process-based calibration metric to do so (Figure 1). We provide optimal gravity drainage parameter values in Table 2. We warn, however, these could be context-specific, as they are known to largely vary depending on simulation context.

640

641

642

643

644

645

2. A high value for the new ice liquid fraction parameter (75%) substantially improves the simulated salinity field (Figure 2, Table 3), especially for ice thinner than 1 m (see Figure 1). This supports the argument that such a value is considered more physically realistic — being more consistent with the absence of salt fractionation during congelation growth (Notz & Worster, 2009), and enabling platelet ice formation (Wongpan et al., 2021).

646

647

648

649

650

3. A vertically-prognostic approach to salt dynamics yields more realistic salinity profiles (Figure 3). This method also displays a more dynamic behavior, in the sense that it increases the variance in vertically-averaged salinity, at all spatial scales — globally, across hemispheres, regionally, and in thickness space — which is more in line with observations.

651

652

653

654

655

656

657

658

4. The type of gravity drainage parameterization remarkably influences the mean salinity, through complex interactions with hydrology, particularly melt pond dynamics (Figure 8). These interaction differ in brine fraction-based (Cox & Weeks, 1988) and Rayleigh-number-based (Griewank & Notz, 2013; Rees Jones & Worster, 2014) parameterizations. By contrast, the two tested Rayleigh number-based parameterizations (GN/RJW) simulate very similar salinity fields, consistent with previous findings in a one-dimensional modeling context (Thomas et al., 2020; Worster & Rees Jones, 2015).
Regarding the impact of vertically-resolved salt dynamics on other sea ice characteristics and on ocean properties, we identify several effects that are broadly consistent with previous studies (Bailey et al., 2020; DuVivier et al., 2021; Turner & Hunke, 2015; Uotila et al., 2017).

Changes in salt dynamics thermodynamically affect the large-scale sea ice mass budget (Figure 7). While inter-simulation differences in sea ice mass budget terms are large, compensation effects mitigate the net sea ice volume response.

1. Using a high value for the new ice liquid fraction parameter $\phi_{new}$ both increases new ice growth in open water (because more saline newly formed ice has reduced enthalpy) and reduces bottom growth (an effect of reduced thermal conductivity).

2. Ice salinity-driven changes in melt pond fraction can induce substantial differences in Arctic sea ice volume between simulations using different gravity drainage parameterizations. This arises from the interplay between the ice salinity profile, surface permeability, pond drainage, albedo, and surface melt rate, as previously identified by Turner & Hunke (2015).

3. Changes in salt dynamics weakly affect the simulated ocean, despite large variations in ice-ocean salt exchanges. This limited oceanic response is caused by compensations between salt drainage and uptake. Also, ice-ocean freshwater fluxes associated with growth and melt invariably accompany salt fluxes and are more influential on the ocean than ice-ocean salt fluxes (Figure 9).

Our analysis leads to several recommendations. First, the optimal value for the new ice liquid fraction parameter is 0.75. Second, we recommend using vertically prognostic salt dynamics combined with a Rayleigh-number-based gravity drainage parameterization, as it gives best overall results in the Arctic (Table 3). Yet, one must remind that in the Antarctic, the diagnostic approach for vertical salinity profile combined with restoring-based gravity drainage parameterization, yields mean salinities and standard deviation that align more closely with observations. This difference may reflect specific, poorly resolved Antarctic salt processes; however, it could also be an artifact of the larger observational uncertainties in the region.

Indeed, the fundamental understanding of sea ice salt dynamics is less advanced in the Antarctic than in the Arctic. In particular, processes important to Antarctic salt dynamics are less well understood and differ markedly from those in the Arctic. These include surface flooding, snow ice formation, and the presence of large near-surface brine tubes (Lewis et al., 2011) in contrast to the brine channels and pockets classically observed in the Arctic (Perovich & Gow, 1996). Moreover, observational salinity biases during core extraction —
potentially more significant in the warmer, more saline Antarctic sea ice — have not been as
thoroughly evaluated as in the Arctic. Finally, our model simulations still largely
underestimate the Antarctic summer ice extent, hence the contribution of multi-year ice to
mean salinity, even in winter, is not accounted for.

In conclusion, NEMO-SI³ (and likely all sea ice models with carefully adjusted salt dynamics) can reproduce observed large-scale ice salinity features. A fully prognostic, vertical treatment improves realism and produces more dynamic salinity variations and ice-ocean exchanges. Salt dynamics strongly affect sea ice mass budget terms, have moderate impact on total sea ice volume, and weakly affect the ocean. Whether they could be important to simulate smaller-scale features (such as fractures, flaw leads, and polynyas) remains to be investigated. Further reducing uncertainties requires deeper understanding of Antarctic sea ice salt dynamics, and better accounting for observational biases.

## Data Availability Statement

The research presented in this paper is based on:

* The NEMO code v5, available freely at https://forge.nemo-ocean.eu/nemo/nemo/-/tree/branch_5.0

* An observational ice core database (BEPSII, 2025), currently unpublished, will also be released on zenodo, regardless of the publication of the paper. This database is also described and discussed in a paper in preparation.

* Five 10-yr-long simulations with NEMO, which monthly, space-dependent output will be released on the zenodo repository, if the paper is published

Meanwhile, the simulations and observations are available on the CNRS cloud, at https://sdrive.cnrs.fr/s/tQRmLnLfb4io6Kb?path=%2F

## Conflict of Interest

The authors declare they have no conflicts of interest relevant to this study.
718 References

719 Abernathey, R. P., Cerovecki, I., Holland, P. R., Newsom, E., Mazloff, M., & Talley, L. D.

720 (2016). Water-mass transformation by sea ice in the upper branch of the Southern

721 Ocean overturning. *Nature Geoscience*, 9(8), 596–601.

722 https://doi.org/10.1038/ngeo2749

723 Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. (2004). Implementation of an Atmosphere–

724 Ocean General Circulation Model on the Expanded Spherical Cube.

725 https://doi.org/10.1175/MWR2823.1

726 Assur, A. (1958). Composition of sea ice and its tensile strength. *A. Nat. Acad. Sci./Nat. Res.*

727 Council, *Arctic Sea Ice*, 598, 106–138.

728 Bailey, D. A., Holland, M. M., DuVivier, A. K., Hunke, E. C., & Turner, A. K. (2020). Impact of

729 a New Sea Ice Thermodynamic Formulation in the CESM2 Sea Ice Component.

730 *Journal of Advances in Modeling Earth Systems*, 12(11), e2020MS002154.

731 https://doi.org/10.1029/2020MS002154

732 Barry, R. G., Serreze, M. C., Maslanik, J. A., & Preller, R. H. (1993). The Arctic Sea Ice–

733 Climate System: Observations and modeling. *Reviews of Geophysics*, 31(4), 397–

734 422. https://doi.org/10.1029/93RG01998

735 BEPSII. (2025). A global compilation of historical temperature and salinity data in sea ice. *To*

736 *Be Released on Zenodo*.

737 Bitz, C. M., & Lipscomb, W. H. (1999). An energy-conserving thermodynamic model of sea

738 ice. *Journal of Geophysical Research*, 104, 15,669–15,677.

739 Bitz, C. M., & Roe, G. H. (2004). A Mechanism for the High Rate of Sea Ice Thinning in the

740 Arctic Ocean. *Journal of Climate*, 17, 3623–3632.

741 Bitz, C. M., Holland, M. M., Weaver, A. J., & Eby, M. (2001). Simulating the ice-thickness

742 distribution in a coupled climate model. *Journal of Geophysical Research*, 106,

743 2441–2463.
744 Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., et al.

745 (2020). Presentation and Evaluation of the IPSL-CM6A-LR Climate Model. *Journal of*

746 *Advances in Modeling Earth Systems*, 12(7), e2019MS002010.

747 https://doi.org/10.1029/2019MS002010

748 Chripko, S., Msadek, R., Sanchez-Gomez, E., Terray, L., Bessières, L., & Moine, M.-P.

749 (2021). Impact of Reduced Arctic Sea Ice on Northern Hemisphere Climate and

750 Weather in Autumn and Winter. https://doi.org/10.1175/JCLI-D-20-0515.1

751 Cox, G. F. N., & Weeks, W. F. (1988). Numerical simulations of the profile properties of

752 undeformed first-year sea ice during the growth season. *Journal of Geophysical*

753 *Research: Oceans*, 93(C10), 12449–12460.

754 https://doi.org/10.1029/JC093iC10p12449

755 Cox, G. F. N., & Weeks, W. S. (1974). Salinity Variations in sea ice. *Journal of Glaciology*,

756 13, 109–120.

757 DuVivier, A. K., Holland, M. M., Landrum, L., Singh, H. A., Bailey, D. A., & Maroon, E. A.

758 (2021). Impacts of Sea Ice Mushy Thermodynamics in the Antarctic on the Coupled

759 Earth System. *Geophysical Research Letters*, 48(18), e2021GL094287.

760 https://doi.org/10.1029/2021GL094287

761 Eicken, H. (1998). Antarctic Sea Ice: Physical Processes, Interactions and Variability. In M.

762 O. Jeffries (Ed.) (Vol. 74, pp. 89–122). American Geophysical Union.

763 Eicken, H., Krouse, H. R., Kadko, D., & Perovich, D. K. (2002). Tracer studies of pathways

764 and rates of meltwater transport through Arctic summer sea ice. *Journal of*

765 *Geophysical Research*, 107. https://doi.org/10.1029/2000JC000583

766 Eicken, Hajo, Gradinger, R., Salganek, M., Shirasawa, K., Perovich, D., & Leppäranta, M.

767 (2010). *Field Techniques for Sea-Ice Research*. *Bibliovault OAI Repository, the*

768 University of Chicago Press.

769 Flocco, D., Feltham, D. L., & Turner, A. K. (2010). Incorporation of a physically based melt

770 pond scheme into the sea ice component of a climate model. *Journal of Geophysical*

771 *Research: Oceans*, 115(C8). https://doi.org/10.1029/2009JC005568
772 Goosse, H., & Fichefet, T. (1999). Importance of ice-ocean interactions for the global ocean circulation: A model study. *Journal of Geophysical Research: Oceans, 104*(C10), 23337-23355. https://doi.org/10.1029/1999JC900215
773
774
775 Griewank, P. J., & Notz, D. (2013). Insights into brine dynamics and sea ice desalination from a 1-D model study of gravity drainage. *Journal of Geophysical Research: Oceans, 118*(7), 3370-3386. https://doi.org/10.1002/jgrc.20247
776
777
778 Holland, M. M., Bailey, D. A., Briegleb, B. P., Light, B., & Hunke, E. (2012). Improved Sea Ice Shortwave Radiation Physics in CCSM4: The Impact of Melt Ponds and Aerosols on Arctic Sea Ice. https://doi.org/10.1175/JCLI-D-11-00078.1
779
780
781 Hunke, E., Allard, R., Bailey, D. A., Blain, P., Craig, A., Dupont, F., et al. (2024, November 19). CICE-Consortium/Icepack: Icepack 1.5.0 (Version 1.5.0). Zenodo. https://doi.org/10.5281/zenodo.14188409
782
783
784 Hunke, E. C., Notz, D. K., Turner, A., & Vancoppenolle, M. (2011). The multi-phase physics of sea ice: a review for modellers. *The Cryosphere*, 5, 989-1009.
785
786 Hunke, Elizabeth C., Hebert, D. A., & Lecomte, O. (2013). Level-ice melt ponds in the Alamos sea ice model, CICE. *Ocean Modelling*, 71, 26-42.
787
788
789 Jeffery, N., Hunke, E., & Elliott, S. (2011). Modeling the transport of passive tracers in sea ice. *Journal of Geophysical Research, under review*.
790
791 Jenkins, M., & Dai, A. (2021). The Impact of Sea-Ice Loss on Arctic Climate Feedbacks and Their Role for Arctic Amplification. *Geophysical Research Letters*, 48(15), e2021GL094599. https://doi.org/10.1029/2021GL094599
792
793
794 Jutras, M., Vancoppenolle, M., Lourenço, A., Vivier, F., Carnat, G., Madec, G., et al. (2016). Thermodynamics of slush and snow-ice formation in the Antarctic sea-ice zone. *Deep Sea Research (II)*, in press.
795
796
797 Kimmritz, M., Danilov, S., & Losch, M. (2016). The adaptive EVP method for solving the sea ice momentum equation. *Ocean Modelling*, 101, 59-67.
798
799
https://doi.org/10.1016/j.ocemod.2016.03.004
800 Kovacs, A. (1996). Sea ice, Part I. Bulk salinity versus ice floe thickness (No. 96–7).

801 Hanover, New Hampshire: Cold Regions Research and Engineering Laboratory.

802 Large, W. G., & Yeager, S. G. (2004). *Diurnal to decadal global forcing for ocean and sea ice models: The datasets and climatologies* (No. TN-460+STR). NCAR Technical Report.

805 de Lavergne, C., Rathore, S., Madec, G., Sallée, J.-B., Ethé, C., Nasser, A., et al. (n.d.).

806 Effects of improved tidal mixing in NEMO one-degree global ocean model. Retrieved from

808 https://www.authorea.com/doi/full/10.22541/essoar.173152139.95978362?commit=5afaf33cf2faa82d5c5d5d43d37fa4f4c5eab210

810 Lavergne, T., Sørensen, A. M., Kern, S., Tonboe, R., Notz, D., Aaboe, S., et al. (2019).

811 Version 2 of the EUMETSAT OSI SAF and ESA CCI sea-ice concentration climate data records. *The Cryosphere*, 13(1), 49–78. https://doi.org/10.5194/tc-13-49-2019

813 Leppäranta, M., Lensu, M., & Witch, P. K. B. (1995). The life story of a first-year sea ice ridge. *Cold Regions Science and Technology*, 23(3), 279–290.

815 Lewis, M. J., Tison, J. L., Weissling, B., Delille, B., Ackley, S. F., Brabant, F., & Xie, H. (2011). Sea ice and snow cover characteristics during the winter-spring transition in the Bellingshausen Sea: An overview of SIMBA 2007. *Deep Sea Research Part II: Topical Studies in Oceanography*, 58(9), 1019–1038.

819 https://doi.org/10.1016/j.dsr2.2010.10.027

820 Light, B., Maykut, G. A., & Grenfell, T. C. (2003). Effects of temperature on the microstructure of first-year Arctic sea ice. *Journal of Geophysical Research*, 108. https://doi.org/10.1029/2001JC000887

823 Lipscomb, W. H. (2001). Remapping the thickness distribution in sea ice models. *Journal of Geophysical Research*, 106, 13989–14000.

825 Madec, G., Bourdallé-Badie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., et al. (2022). NEMO ocean engine. https://doi.org/10.5281/zenodo.6334656
827 Malmgren, F. (1927). The Norwegian North Polar Expedition with the "Maud" 1918-1925 vol.
828 1a no.5. In H. U. Sverdrup (Ed.) (pp. 1-67). John Griegds Boktr., Bergen, Norway.
829 Manabe, S., & Stouffer, R. J. (1980). Sensitivity of a Global Climate Model to an Increase of
830 CO₂ Concentration in the Atmosphere. *Journal of Geophysical Research*, 85, 5529-
831 5554.
832 Notz, D. (2005). *Thermodynamic and fluid-dynamical processes in sea ice*. University of
833 Cambridge.
834 Notz, D., & Worster, M. G. (2008). In situ measurements of the evolution of young sea ice.
835 *Journal of Geophysical Research: Oceans*, 113(C3).
836 https://doi.org/10.1029/2007JC004333
837 Notz, D., & Worster, M. G. (2009). Desalination processes of sea ice revisited. *Journal of
838 Geophysical Research: Oceans*, 114(C5). https://doi.org/10.1029/2008JC004885
839 Notz, D., Wettlaufer, J. S., & Worster, M. G. (2005). A non-destructive method for measuring
840 the salinity and solid fraction of growing sea ice in situ. *Journal of Glaciology*, 51,
841 159-166.
842 Ono, N. (1967). Specific heat and heat of fusion of sea ice. In H. Oura (Ed.), *Physics of snow
843 and ice: Proceedings of the International Conference on Low Temperature Science*.
844 *I. Conference on physics of snow and ice*, (pp. 599-610). Institute of Low
845 Temperature Science, Hokkaido, Japan.
846 Peixoto, J. P., & Oort, A. H. (1992). *Physics of climate*. (Springer-Verlag, Ed.). Springer-
847 Verlag.
848 Pellichero, V., Sallée, J.-B., Chapman, C. C., & Downes, S. M. (2018). The southern ocean
849 meridional overturning in the sea-ice sector is driven by freshwater fluxes. *Nature
850 Communications*, 9(1), 1789. https://doi.org/10.1038/s41467-018-04101-2
851 Perovich, D. K., & Gow, A. J. (1996). A quantitative description of sea ice inclusions. *Jounal
852 of Geophysical Research*, 101, 18327-18343.
853 Prather, M. J. (1986). Numerical advection by conservation of second-order moments.
854 *Journal of Geophysical Research*, 91, 6671-6681.
855 Pringle, D. J., Eicken, H., Trodahl, H. J., & Backstrom, L. G. E. (2007). Thermal conductivity of landfast Antarctic and Arctic sea ice. *Journal of Geophysical Research*, 112, doi:10.1029/2006JC003641. https://doi.org/10.1029/2006JC003641

856
857

858 Rees Jones, D. W., & Worster, M. G. (2014). A physically based parameterization of gravity drainage for sea-ice modeling. *Journal of Geophysical Research: Oceans*, 119(9), 5599–5621. https://doi.org/10.1002/2013JC009296

859
860

861 Roach, L. A., Dörr, J., Holmes, C. R., Massonnet, F., Blockley, E. W., Notz, D., et al. (2020). Antarctic Sea Ice Area in CMIP6. *Geophysical Research Letters*, 47(9), e2019GL086729. https://doi.org/10.1029/2019GL086729

862
863

864 Schmidt, G., Bitz, C. M., Mikolajewicz, U., & Tremblay, L.-B. (2004). Ice-ocean boundary conditions for coupled models. *Ocean Modelling*, 7, 59–74.

865
866

867 SIMIP Community. (2020). Arctic Sea Ice in CMIP6. *Geophysical Research Letters*, 47(10), e2019GL086749. https://doi.org/10.1029/2019GL086749

868
869

870 Smedsrud, L. H., & Skogseth, R. (2006). Field measurements of Arctic grease ice properties and processes. *Cold Regions Science and Technology*, 44(3), 171–183.

871
872

873 Sterlin, J., Fichefet, T., Massonnet, F., Lecomte, O., & Vancoppenolle, M. (2021). Sensitivity of Arctic sea ice to melt pond processes and atmospheric forcing: A model study. *Ocean Modelling*, 167, 101872. https://doi.org/10.1016/j.ocemod.2021.101872

874
875

876 Stössel, A., Yang, K., & Kim, S. J. (2002). On the role of sea ice and convection in a global ocean model. *Journal of Physical Oceanography*, 32, 1194–1209.

877
878

879 Stroeve, J., & Notz, D. (2018). Changing state of Arctic sea ice across all seasons. *Environmental Research Letters*, 13(10), 103001. https://doi.org/10.1088/1748-9326/aade56

880
881

882 Thomas, M., Vancoppenolle, M., France, J., Sturges, W., Bakker, D., Kaiser, J., & von Glasow, R. (2020). Tracer measurements in growing sea ice support convective gravity drainage parameterisations. *Journal of Geophysical Research*.
883 Thorndike, A. S., Rothrock, D. A., Maykut, G. A., & Colony, R. (1975). The thickness distribution of sea ice. *Journal of Geophysical Research*, 80, 4501–4513.

884 Toyoda, T., Sakamoto, K., Toyota, T., Tsujino, H., Urakawa, L. S., Kawakami, Y., et al. (2024). Improvement of sea ice thermodynamics with variable sea ice salinity and melt pond parameterizations in an OGCM. *Ocean Modelling*, 187, 102288.
888 https://doi.org/10.1016/j.ocemod.2023.102288

889 Turner, A. K., & Hunke, E. C. (2015). Impacts of a mushy-layer thermodynamic approach in global sea-ice simulations using the CICE sea-ice model. *Journal of Geophysical Research: Oceans*, 120(2), 1253–1275. https://doi.org/10.1002/2014JC010358

892 Turner, A. K., Hunke, E. C., & Bitz, C. M. (2013). Two modes of sea-ice gravity drainage: A parameterization for large-scale modeling. *Journal of Geophysical Research: Oceans*, 118(5), 2279–2294. https://doi.org/10.1002/jgrc.20171

895 Untersteiner, N. (1968). Natural desalination and equilibrium salinity profile of perennial sea ice. *Journal of Geophysical Research*, 73, 1251–1257.

897 Uotila, P., Iovino, D., Vancoppenolle, M., Lensü, M., & Rousset, C. (2017). On the influence of sea-ice physics in multi-decadal ocean-ice hindcasts. *Geoscientific Model Development*, 10–2, 1009–1031. https://doi.org/10.5194/gmd-2016-187

900 Vancoppenolle, M. (2023). SI3, the NEMO Sea Ice Engine. https://doi.org/10.5281/zenodo.7534900

902 Vancoppenolle, M., Fichefet, T., & Bitz, C. M. (2005). On the sensitivity of undeformed Arctic sea ice to its vertical salinity profile. *Geophysical Research Letters*, 32. https://doi.org/10.1029/2005GL023427

905 Vancoppenolle, M., Fichefet, T., & Bitz, C. M. (2006). Modeling the salinity profile of undeformed Arctic sea ice. *Geophysical Research Letters*, 33. https://doi.org/10.1029/2006GL028342

908 Vancoppenolle, M., Fichefet, T., & Goosse, H. (2009). Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 2. Sensitivity to salinity processes. *Ocean Modelling*, 27 (1–2), 54–69. https://doi.org/10.1016/j.ocemod.2008.11.003
911 Vancoppenolle, Martin, Bitz, C. M., & Fichefet, T. (2007). Summer landfast sea ice

912 desalination at Point Barrow: Model and Observations. *Journal of Geophysical*

913 *Research*, 112. https://doi.org/10.1029/2006JC003493

914 Vancoppenolle, Martin, Goosse, H., Montety, A. de, Fichefet, T., Tremblay, B., & Tison, J.-L.

915 (2010). Modeling brine and nutrient dynamics in Antarctic sea ice: the case of

916 dissolved silica. *Journal of Geophysical Research*, 115 (C2).

917 https://doi.org/10.1029/2009JC005369

918 Vancoppenolle, Martin, Madec, G., Thomas, M., & McDougall, T. J. (2019).

919 Thermodynamics of sea ice phase composition revisited. *Journal of Geophysical*

920 *Research: Oceans*, 0(ja). https://doi.org/10.1029/2018JC014611

921 Weeks, W. F., & Ackley, S. F. (1986). The Geophysics of sea ice. In N. Untersteiner (Ed.)

922 (Vol. 146, pp. 9–164). Plenum, New York.

923 Wells, A. J., Wettlaufer, J. S., & Orszag, S. A. (2011). Brine fluxes from growing sea ice.

924 *Geophys. Res. Lett.*, 38. https://doi.org/10.1029/2010GL046288

925 Wongpan, P., Vancoppenolle, M., Langhorne, P. J., Smith, I. J., Madec, G., Gough, A. J., et

926 al. (2021). Sub-Ice Platelet Layer Physics: Insights From a Mushy-Layer Sea Ice

927 Model. *Journal of Geophysical Research: Oceans*, 126(6), e2019JC015918.

928 https://doi.org/10.1029/2019JC015918

929 Worster, M. G., & Rees Jones, D. W. (2015). Sea-ice thermodynamics and brine drainage.

930 *Phil. Trans. R. Soc. A*, 373 (20140166). https://doi.org/10.1098/rsta.2014.0166

931
<table><thead><tr><td>Simulations</td><td>Approach</td><td>Main formula</td><td>Reference</td></tr></thead><tbody><tr><td>RESTO_REF, RESTO_MOD</td><td>(1) Restoring</td><td>d<i>S</i>(t)/dt = -(<i>S</i> - <i>S</i><sub>∞</sub>)/<i>τ</i></td><td>Vancoppenolle et al (2009)</td></tr><tr><td>CW</td><td>(2) Ramping function of brine fraction</td><td><i>w</i><sub>CW</sub>(<i>z</i>) = -<i>μ</i> &middot; <i>α</i><sub>CW</sub> &middot; max(<sup><i>φ</i>(<i>z</i>)</sup>&frasl;<sub><i>φ</i><sub>c</sub></sub> - 1, 0)</td><td>Cox and Weeks (1988), Thomas et al (2020)</td></tr><tr><td>GN</td><td>(3) Rayleigh number functional (integral)</td><td><i>w</i><sub>GN,k</sub> = -<sup><i>α</i><sub>GN</sub></sup>&frasl;<sub><i>ρ</i></sub> &sum;<sub><i>i</i>=1</sub><sup><i>k</i></sup> (<i>Ra</i><sub><i>i</i></sub> - <i>Ra</i><sub>c</sub>)&Delta;<i>z</i><sub><i>i</i></sub></td><td>Griewank and Notz (2013)</td></tr><tr><td>RJW</td><td>(4) Rayleigh number functional (max)</td><td><i>w</i><sub>RJW</sub>(<i>z</i>) = -<i>α</i><sub>RJW</sub> &middot; <sup><i>k</i></sup>&frasl;<sub><i>c</i></sub> max<sub><i>z</i><sub>c</sub> &le; <i>z</i> &le; <i>h</i></sub> [<i>Ra</i>(<i>z</i>) - <i>Ra</i><sub>c</sub>] &middot; <sup><i>z</i> - <i>z</i><sub>c</sub></sup>&frasl;<sub>(<i>h</i> - <i>z</i><sub>c</sub>)<sup>2</sup></sub></td><td>Rees Jones and Worster (2014)</td></tr></tbody></table>

933

**Table 1.** Details of the gravity drainage parameterizations used in this work. Symbols are defined as follows. $\hat{S}$ is the vertically-averaged ice salinity, $w$ is the convective velocity, $\phi$ is brine fraction, $Ra$ is the Rayleigh number, and $Ra_c$ its critical value. $\alpha_{CW}$ [g·m·(kg·°C·s)<sup>-1</sup>] $\alpha_{GN}$ [kg·(m3·s)<sup>-1</sup>] and $\alpha_{RJW}$ (dimensionless) are the respective brine convection intensity tuning parameters in each of the parameterizations. In parameterization (1), $S_\infty = 5$ g/kg is the asymptotic cold thick ice salinity (g/kg) and $\tau$ is the gravity drainage restoring time scale (s). In (2), $\mu = -0.054$ °C·kg·g<sup>-1</sup> is the freezing temperature depression coefficient and $\phi_c = 0.05$ is the critical porosity threshold for permeability. Convective velocity is non-zero only if the temperature gradient points upwards. In (3), $\rho$ is brine density and $\Delta z_i$ is the thickness of the ith vertical layer. In (4) $c$, $k$ are the brine heat capacity and thermal conductivity and $z_c$ is the critical depth, i.e. the uppermost depth in the sea ice where $Ra$ is above $Ra_c$. Parameterizations (3) and (4) both assume a zero Darcy velocity above $z = z_c$.

946
<table><caption>Table 2. Parameterization choices for sea ice salt dynamics in the simulations used in this work: value of the new ice liquid fraction parameter ($\phi_{new}$); choice of a diagnostic or prognostic vertical salinity profile; gravity drainage parameterization retained and parameter values. For comparison, parameter values used in Thomas et al (2020, T20), set in the context of artificial young sea ice simulations, are also given.</caption><thead><tr><th rowspan="2">Simulation ID</th><th rowspan="2">φ<sub>new</sub></th><th rowspan="2">Approach for vertical profile</th><th colspan="4">Gravity drainage parameterization</th></tr><tr><th>Approach and reference</th><th>Parameter(s) symbol and units</th><th>Parameter values (this study)</th><th>Parameter values (T20)</th></tr></thead><tbody><tr><td>RESTO_REF</td><td>0.25</td><td>Diagnostic</td><td>Salinity restoring (Vancoppenolle et al, 2009)</td><td>τ [days]</td><td>20</td><td>n.a.</td></tr><tr><td>RESTO_MOD</td><td>0.75</td><td>Diagnostic</td><td>Salinity restoring (Vancoppenolle et al, 2009)</td><td>τ [days]</td><td>2</td><td>n.a.</td></tr><tr><td>CW</td><td>0.75</td><td>Prognostic</td><td>Liquid fraction-based Darcy velocity (Cox and Weeks, 1988)</td><td>α<sub>CW</sub> [g·m·(kg·°C·s)<sup>-1</sup>]</td><td>4.8 x 10<sup>-7</sup></td><td>7.2 x 10<sup>-7</sup></td></tr><tr><td>GN</td><td>0.75</td><td>Prognostic</td><td>Ra-number-based Darcy velocity (Griewank and Notz, 2013)</td><td>Ra<sub>c</sub>, α<sub>GN</sub> [kg·(m<sup>3</sup>·s)<sup>-1</sup>]</td><td>10, 5.1 x 10<sup>-4</sup></td><td>2.4, 6.7 x 10<sup>-3</sup></td></tr><tr><td>RJW</td><td>0.75</td><td>Prognostic</td><td>Ra-number based Darcy velocity (Rees Jones and Worster, 2015)</td><td>Ra<sub>c</sub>, α<sub>RJW</sub> [.]</td><td>9 0.03</td><td>2.9 0.13</td></tr></tbody></table>
<table><caption>**Table 3.** Hemispheric salinity statistics across the different simulations: mean salinity (<S>), standard deviation (<σ_S>) and low-surface salinity probability (f<sub>low-S</sub>, i.e. the fraction of cores / grid cells with a salinity of the upper 20 cm being less than 2 g/kg). Source data for simulations are monthly mean fields of vertically-averaged salinity, S(x,y,month), for all grid cells with sea ice concentration > 15%. Source observational data include all available ice cores from the corresponding hemisphere and sampling month.</caption><thead><tr><th rowspan="2">Simulation</th><th colspan="3">Northern hemisphere</th><th colspan="3">Southern hemisphere</th><th rowspan="2">SH-NH<br>Δ&lt;S&gt;</th></tr><tr><th>&lt;S&gt;<br>(g/kg)</th><th>&lt;σ<sub>S</sub>&gt;<br>(g/kg)</th><th>f<sub>low-S</sub><br>(%)</th><th>&lt;S&gt;<br>(g/kg)</th><th>&lt;σ<sub>S</sub>&gt;<br>(g/kg)</th><th>f<sub>low-S</sub><br>(%)</th></tr></thead><tbody><tr><td>RESTO_REF</td><td>3.8</td><td>0.9</td><td>45</td><td>5.3</td><td>1.0</td><td>6.8</td><td>1.5</td></tr><tr><td>RESTO_MOD</td><td>4.4</td><td>1.4</td><td>35</td><td>6.2</td><td>2.8</td><td>6.9</td><td>1.9</td></tr><tr><td>CW</td><td>5.1</td><td>1.7</td><td>51</td><td>8.0</td><td>2.9</td><td>3.1</td><td>2.9</td></tr><tr><td>GN</td><td>4.3</td><td>1.7</td><td>54</td><td>7.6</td><td>2.9</td><td>3.1</td><td>3.3</td></tr><tr><td>RJW</td><td>4.4</td><td>1.7</td><td>49</td><td>7.8</td><td>2.9</td><td>1.8</td><td>3.4</td></tr><tr><td>OBS</td><td>4.4</td><td>1.8</td><td>35</td><td>6.6</td><td>2.7</td><td>7.3</td><td>2.2</td></tr></tbody></table>
<figure><img src="image_3.png" /></figure>

960

961 **Figure 1.** Geographical distribution of core extraction location. Colors indicate sampling

962 decade and symbol type the sampling season. Shading shows the passive microwave

963 satellite 1979-2020 average sea ice concentration, at the month of ice extent maximum —

964 March in the Arctic, September in the Antarctic (Lavergne et al., 2019).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Ice salinity-thickness relationship in the five simulations (**a-e**; year 10 of the simulation), and in the observational data (**f**, dots). Grid-cells / cores with upper-20 cm ice salinity less than 2 g/kg were excluded from analysis. Model grid cells with less than 50% ice concentration were also excluded. Lines are $y = a + b/x$ fits, based on either observations (black lines) or simulations (colored lines). Colored shading in panels **a-e** displays 2D density histograms of simulated mean ice salinity (inter-category mean value, weighted by ice fraction). The gray shading in panels **a-f** indicates the envelope defined by the observed minimum and maximum values, within log-uniform thickness bins.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Simulated versus observed seasonal cycle of monthly mean sea ice salinity in the (a) Northern and (b) Southern hemispheres. Solid lines give simulated ice-concentration-weighted mean salinity. Grid cells with ice concentration under 50% were excluded from the means. Box plots give corresponding observed monthly mean and statistics (box: 10-90 percentiles; black horizontal line: mean). The number of observations for each month is given in the grey box. Fall, winter, spring, and summer seasons are respectively defined as Oct-Dec, Jan-Mar, Apr-Jun, and Jul-Sep for the Arctic; and Mar-May, Jun-Aug, Sep-Nov, and Dec-Feb for the Antarctic.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mean depth-normalized vertical sea ice salinity profile for each hemisphere (top: Northern hemisphere; bottom: Southern hemisphere) and within each surface-salinity class, i.e. high-surface salinity (HSS, left) and low-surface salinity (LSS, right) classes. Observations are in black, with the 10th to 90th percentiles shaded in gray. Colored solid lines give simulated mean profiles. The number of cores used in the observed mean profile calculation is given. Cores with less than three sections were excluded. Model data include only grid cells where sea ice concentration exceeds 50%, and means are weighted by sea ice area fraction)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Geographical distribution of sea ice salinity, along with sea ice thickness contours (green is 1 m, purple 2 m and black 3 m) in the different simulations (RESTO_REF, RESTO_MOD, CW, GN, and RJW) for the Northern hemisphere in March (a-d) and September (e-h) 2009, and for the Southern hemisphere in September (i-l) 2009. Summer maps for the Southern hemisphere are not shown because there is very little ice left in our simulations.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Simulated seasonal cycles of sea ice volume in the (a) Northern and (b) Southern hemisphere, for the 10th year of the five simulations analyzed in this paper. Insets give annual mean volume versus annual mean salinity.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

1005

1006 **Figure 7.** Annually and spatially integrated ocean-to-ice mass and salt fluxes over the (a, b)

1007 Northern and (c, d) Southern hemispheres, diagnosed from year 10 of the different
1008 simulations. The mass budget is decomposed into all growth and melt processes, whereas
1009 the salt budget is split into the contributions of brine drainage, ice growth and melt (i.e.
1010 growth and melt sub-processes are grouped). Positive values indicate a gain for sea ice.

1011
<figure><img src="image_13.png" /></figure>

1012

1013 **Figure 8.** Analysis of the causes of sea ice thickness differences in the CW simulation as compared with GN or RJW simulations, during the first simulation year. **a)** Seasonal cycles of the Northern hemisphere simulated sea ice volume in four experiments (RESTO_MOD, CW, GN and RJW). **b-f)** Spatial distribution of the CW minus GN differences, averaged over June-July of the first simulation year: **b)** vertically-averaged ice salinity, **c)** vertically-averaged ice temperature, **d)** mean sea ice thickness, **e)** surface melt rate, **f)** melt pond fraction. Analysis is similar if RJW is used instead of GN.

1020
<figure><img src="image_14.png" /></figure>

1021

1022 **Figure 9.** Net ocean-to-ice salt flux (top), and estimated salinity tendency of the upper ocean (bottom), both plotted against water volume fluxes in the different simulations, based on monthly mean gridded outputs, for each month of year 10. The salinity tendency of the upper ocean is estimated for a 20 m surface layer with a salinity of 34 g/kg. NH (SH) data are shown in red (blue).