# Tracer measurements in growing sea ice support convective gravity drainage parameterizations

M. Thomas, Martin Vancoppenolle, J.-L. France, W. T. Sturges, Dorothee C.

E. Bakker, J. Kaiser, R. von Glasow

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

## To cite this version:

M. Thomas, Martin Vancoppenolle, J.-L. France, W. T. Sturges, Dorothee C. E. Bakker, et al., Tracer measurements in growing sea ice support convective gravity drainage parameterizations. Journal of Geophysical Research. Oceans, 2020, 125 (2), pp.e2019JC015791. 10.1029/2019jc015791. hal-02904132

HAL Id: hal-02904132

https://hal.science/hal-02904132

Submitted on 30 Nov 2020

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>
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

10.1029/2019JC015791

**Key Points:**

* We use measurements of bulk salinity and rhodamine concentrations to trace brine dynamics in a sea ice growth experiment

* Convective gravity drainage parameterizations using Rayleigh numbers perform best against observational metrics

* These convective parameterizations are appropriate for predicting the dynamics of any solute in growing sea ice

**Supporting Information:**

* Figure S1

* Figure S2

* Figure S3

**Correspondence to:**

M. Thomas,
max.thomas@uea.ac.uk

**Citation:**

Thomas, M., Vancoppenolle, M., France, J. L., Sturges, W. T., Bakker, D. C. E., Kaiser, J., & von Glasow, R. (2020). Tracer measurements in growing sea ice support convective gravity drainage parameterizations. *Journal of Geophysical Research: Oceans*, 125, e2019JC015791. https://doi.org/10.1029/2019JC015791

Received 23 OCT 2019

Accepted 13 JAN 2020

Accepted article online 22 JAN 2020

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

# Tracer Measurements in Growing Sea Ice Support Convective Gravity Drainage Parameterizations

M. Thomas¹ ID, M. Vancoppenolle² ID, J. L. France¹,³,⁴, W. T. Sturges¹, D. C. E. Bakker¹ ID, J. Kaiser¹ ID, and R. von Glasow¹

¹Centre for Ocean and Atmospheric Sciences, School of Environmental Sciences, University of East Anglia, Norwich, UK, ²Sorbonne Université, Laboratoire d'Océanographie et du Climat (LOCEAN), Institut Pierre-Simon-Laplace (IPSL), CNRS/IRD/MNHN, Paris, France, ³British Antarctic Survey, Cambridge, UK, ⁴Department of Earth Sciences, Royal Holloway, University of London, Egham, UK

**Abstract** Gravity drainage is the dominant process redistributing solutes in growing sea ice. Modeling gravity drainage is therefore necessary to predict physical and biogeochemical variables in sea ice. We evaluate seven gravity drainage parameterizations, spanning the range of approaches in the literature, using tracer measurements in a sea ice growth experiment. Artificial sea ice is grown to around 17 cm thickness in a new experimental facility, the Roland von Glasow air-sea-ice chamber. We use NaCl (present in the water initially) and rhodamine (injected into the water after 10 cm of sea ice growth) as independent tracers of brine dynamics. We measure vertical profiles of bulk salinity in situ, as well as bulk salinity and rhodamine in discrete samples taken at the end of the experiment. Convective parameterizations that diagnose gravity drainage using Rayleigh numbers outperform a simpler convective parameterization and diffusive parameterizations when compared to observations. This study is the first to numerically model solutes decoupled from salinity using convective gravity drainage parameterizations. Our results show that (1) convective, Rayleigh number-based parameterizations are our most accurate and precise tool for predicting sea ice bulk salinity; and (2) these parameterizations can be generalized to brine dynamics parameterizations, and hence can predict the dynamics of any solute in growing sea ice.

**Plain Language Summary** The cold atmosphere in the Arctic and Southern Oceans can cause sea water to freeze at the surface, forming sea ice. The salt present in the sea water is trapped in the newly formed sea ice, in pockets of very salty liquid, or brine. Brine is much denser than the ocean below, which causes the brine to drain into the ocean. This process is generally referred to as gravity drainage. It is critically important, affecting the ocean's circulation, the movement of greenhouse gases through sea ice, and the supply of nutrients to sea ice algae. Several authors have suggested ways to model gravity drainage. Broadly, these can be split into two approaches: convective and diffusive. Here, we use a laboratory sea ice growth experiment to test these models. We track the movement of brine using measurements of two chemicals in the sea ice and compare measured concentrations to those predicted by the models. Our results show that convective approaches are our most powerful tool for predicting the movement of brine by gravity drainage and do an excellent job.

## 1. Introduction

Vast areas of the world's oceans freeze each year. Some of the salt originally contained in the freezing sea-water is retained in the resulting sea ice, concentrated into liquid brine inclusions, and some is lost to the ocean below. The dominant process redistributing brine within growing sea ice and exchanging brine with the ocean is gravity drainage (Notz & Worster, 2009). As such, gravity drainage exerts a control on the habitable volume (Vancoppenolle et al., 2013), nutrient budget (Fripiat et al., 2017; Fritsen et al., 1994), and inorganic carbon budget (König et al., 2018; Moreau et al., 2015) in growing sea ice, as well as determining the buoyancy forcing provided by growing sea ice to the ocean (Pellichero et al., 2018; Worster & Rees Jones, 2015).

Gravity drainage is driven by the greater density of brine relative to the ocean below. In growing sea ice the coldest (and therefore most saline and densest) brine lies above less dense brine and ocean (Cox & Weeks, 1975). Laboratory (Eide & Martin, 1975; Middleton et al., 2016; Niedrauer & Martin, 1979) and field (Lake

<footnote>©2020. The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits
use, distribution and reproduction in
any medium, provided the original
work is properly cited.</footnote>
& Lewis, 1970) observations have identified a downward flow of brine through narrow, vertically aligned drainage areas called “brine channels.” The liquid fraction of the bulk sea ice moderates the flow of brine by setting the sea ice permeability (Freitag, 1999). The earliest, simplest parameterization of gravity drainage (Cox & Weeks, 1988) empirically relates the vertical temperature gradient and liquid fraction to the sea ice/ocean salt flux.

Notz and Worster (2008) documented the evolution of sea ice salinity during the early stages of natural Arctic sea ice growth using nondestructive salinity measurements (Notz et al., 2005). Their analysis suggests that gravity drainage intensity relates to a depth-dependent Rayleigh number describing the competition between buoyant instability and dissipation (Rees Jones & Worster, 2014). High-resolution two-dimensional numerical simulations of growing sea ice based on mushy-layer equations (Wells et al., 2011) indicate a linear relationship between the salt flux from growing sea ice due to gravity drainage and the Rayleigh number. These simulations, and earlier experimental evidence (Eide & Martin, 1975), also suggest the existence of an upward return flow of brine over a broad mush region, replacing the downward flow in narrow brine channels. In this context, gravity drainage can be seen as the convective overturning of brine.

Following these works, several one-dimensional, vertically resolved parameterizations were proposed seeking to improve on Cox and Weeks (1988) and possible inclusion into Earth System Models. Griewank and Notz (2013), Turner et al. (2013), and Rees Jones and Worster (2014) parameterize gravity drainage as a convective process and formulate the velocity of moving brine using the Rayleigh number. Worster and Rees Jones (2015) show that these parameterizations are broadly similar. Vancoppenolle et al. (2010) parameterize gravity drainage as a diffusive process, using a Rayleigh number to set the strength of the diffusion. Jeffery et al. (2011) propose two diffusive parameterizations. One uses a mixing length, similar in form to the Rayleigh number, to set the strength of diffusion, while the other uses only the liquid fraction. No quantitative comparison of these salt dynamics parameterizations has yet been performed.

Another motivation for the development of improved gravity drainage parameterizations was the need to simulate solutes other than salt. Indeed, the classical, empirical approach of Cox and Weeks (1988) directly predicts desalination rate and cannot readily be applied to solutes in general. A similar framework for both salt and other solutes (Jeffery et al., 2011; Vancoppenolle et al., 2010) was applied to study the mass budget of biogeochemical compounds in sea ice (Kotovitch et al., 2016; Moreau et al., 2014, 2015). To date, though physical arguments favor convective parameterizations (Worster & Rees Jones, 2015), only diffusive schemes have been used to represent the transport of solutes that are decoupled from the salinity by gravity drainage.

The transport of solutes in general is a priori governed by the same processes that transport salt. However, it does not necessarily follow that gravity drainage parameterizations, tuned and evaluated for salt, would work as well for other solutes. One major specificity of salt is that brine salinity is practically determined by temperature (Assur, 1958; Worster, 1992) because of the quick adjustment toward thermal equilibrium (Feltham et al., 2006; Jutras et al., 2016; Worster, 1992). By contrast, many other solutes can be created or destroyed in sea ice by biological or chemical processes, decoupling their brine concentrations from brine salinity (see, e.g., Fripiat et al., 2017, for the case of nutrients). Whereas the vertical brine salinity profile in growing sea ice is monotonic, such a statement cannot be made for solutes in general. These arguments suggest the need to verify that gravity drainage parameterizations, supposedly appropriate to salt dynamics, are suitable for any solute. Doing so would allow us to confidently generalize these salt dynamics parameterizations to brine dynamics parameterizations.

Here, we address the following questions: (1) How well do current gravity drainage parameterizations predict salt dynamics in sea ice and which are best; and (2) How well do these parameterizations predict brine dynamics, and hence the dynamics of any solute?

We do so by comparing the output of seven gravity drainage parameterizations, covering the range of approaches in the literature, to measurements of two independent tracers of brine dynamics. The tracers used were NaCl, and rhodamine 6G which was deliberately decoupled from the salinity. Artificial sea ice was grown in a new experimental facility, the Roland von Glasow air-sea-ice chamber.

## 2. Modeling Gravity Drainage as a 1-D Vertical Process

In this section, we briefly describe the one-dimensional vertical gravity drainage parameterizations that have been proposed. All simulate the vertical bulk salinity profile as a function of time, $S(z, t)$. For all
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic diagram of the transport of solutes in (left) convective and (right) diffusive gravity drainage parameterizations. In the convective scheme (adapted from Griewank & Notz, 2013), parcels of brine are moved from each layer directly into the ocean, and this brine is replaced by the adjacent underlying layer, or the ocean. In the diffusive scheme solutes diffuse between sea ice layers, and between the sea ice and the ocean.</figcaption></figure>

of them, thermodynamic equilibrium is assumed, such that brine salinity is constrained by temperature ($S_{br} = S_{br}(T)$).

## 2.1. Convection Diagnosed by the Rayleigh Number

Several convective one-dimensional gravity drainage parameterizations have been proposed (Griewank & Notz, 2013; Rees Jones & Worster, 2014; Turner et al., 2013). Convective schemes derive from ideas developed in Wells et al. (2011) with a downward brine channel flow from each layer directly connecting the sea ice interior to the ocean, and a compensating upward return flow (Figure 1). These schemes essentially solve a one-dimensional salt advection equation of the form

$$ \frac{\partial S}{\partial t} = -w \frac{\partial S_{br}}{\partial z}, \qquad (1) $$

where $w$ is the upwelling Darcy velocity of the return flow (positive downward; Rees Jones & Worster, 2014), $S_{br}$ is the brine salinity, and $z$ and $t$ are depth and time, respectively. Both $w$ and the brine salinity gradient are negative for growing sea ice, giving a net desalination. The form of equation (1), intentionally not conservative as the salt lost from each sea ice layer, is assumed to be deposited directly into the underlying water via brine channels.

Rees Jones and Worster (2014), from here referred to as “RJW14,” build on previous analysis (Rees Jones & Worster, 2013a, 2013b), to formulate $w$ as follows. First, they assume that brine convection occurs between a critical depth, $z_c$, and the ice-ocean interface. $z_c$ is defined as the uppermost depth in the sea ice where the Rayleigh number, $Ra$, reaches some critical value, $Ra_c$. From Rees Jones and Worster (2014), $Ra$ is given by

$$ Ra(z) = \frac{cg\beta(S_{br}(z) - S_w)\bar{\Pi}}{kv} (h - z), \qquad (2) $$

where $z$ is depth counted from the upper sea ice surface and is positive downward and $h$ is sea ice thickness. $S_w$ is the ocean salinity and $\beta$ is the saline density coefficient. $\bar{\Pi}$ is the effective permeability, which is a priori depth dependent. $g$ is the acceleration due to gravity, and $c$, $k$, and $v$ are the heat capacity, thermal conductivity, and kinematic viscosity of brine, respectively (Table 1). Second, RJW14 assume that the intensity of convection is proportional to the maximum Rayleigh number over the convective layer. Within the convecting layer, RJW14 use a linear dependence of $w$ with depth,

$$ w = -\alpha_{RJW} \cdot \frac{k}{c} \max_{z_c \le z \le h} [Ra(z) - Ra_c] \cdot \frac{z - z_c}{(h - z_c)^2}. \qquad (3) $$

$\alpha_{RJW}$ is a dimensionless convection strength tuning parameter and $\rho$ is the density of brine. Outside of the convecting layer, $w$ is zero.
<table><caption>**Table 1**
Symbols Used in This Work</caption><thead><tr><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>A</td><td>internal sea ice area density</td><td>cm<sup>2</sup>/kg</td></tr><tr><td>c</td><td>brine heat capacity</td><td>J·m<sup>-3</sup>·K<sup>-1</sup></td></tr><tr><td>C</td><td>bulk rhodamine concentration</td><td>µg/kg</td></tr><tr><td>C<sub>ads</sub></td><td>rhodamine adsorbed to<br>internal surfaces</td><td>µg/kg</td></tr><tr><td>C<sub>br</sub></td><td>brine rhodamine concentration</td><td>µg/kg</td></tr><tr><td>C<sub>w</sub></td><td>underlying water rhodamine<br>concentration</td><td>µg/kg</td></tr><tr><td>D<sub>e</sub></td><td>effective diffusivity</td><td>m<sup>2</sup>/s</td></tr><tr><td>D<sub>tur</sub></td><td>turbulent diffusivity</td><td>m<sup>2</sup>/s</td></tr><tr><td>D<sub>m</sub></td><td>molecular diffusivity</td><td>m<sup>2</sup>/s</td></tr><tr><td>F<sub>C</sub></td><td>rhodamine flux at upper interface</td><td>µg·m<sup>-2</sup>·s<sup>-1</sup></td></tr><tr><td>F<sub>S</sub></td><td>salt flux at upper interface</td><td>g·m<sup>-2</sup>·s<sup>-1</sup></td></tr><tr><td>g</td><td>acceleration due to gravity</td><td>m/s<sup>2</sup></td></tr><tr><td>h</td><td>sea ice thickness</td><td>m</td></tr><tr><td>k</td><td>brine thermal conductivity</td><td>W·m<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>l<sub>m</sub></td><td>mixing length</td><td>m</td></tr><tr><td>m<sub>br↑</sub></td><td>upwelling brine mass</td><td>kg/m<sup>2</sup></td></tr><tr><td>m<sub>br↓</sub></td><td>downwelling brine mass</td><td>kg/m<sup>2</sup></td></tr><tr><td>Ra</td><td>Rayleigh number</td><td></td></tr><tr><td>Ra<sub>c</sub></td><td>critical Rayleigh number</td><td></td></tr><tr><td>S</td><td>bulk sea ice salinity</td><td>g/kg</td></tr><tr><td>S<sub>br</sub></td><td>brine salinity</td><td>g/kg</td></tr><tr><td>S<sub>w</sub></td><td>underlying water salinity</td><td>g/kg</td></tr><tr><td>t</td><td>time</td><td>s</td></tr><tr><td>T</td><td>sea ice temperature</td><td>°C</td></tr><tr><td>T<sub>s</sub></td><td>temperature at upper interface</td><td>°C</td></tr><tr><td>w</td><td>upwelling brine velocity</td><td>m/s</td></tr><tr><td>z</td><td>depth in sea ice</td><td>m</td></tr><tr><td>z<sub>c</sub></td><td>depth of top of convecting layer</td><td>m</td></tr><tr><td>Z</td><td>resistance between wirepair</td><td>Ω</td></tr><tr><td>Z<sub>0</sub></td><td>resistance between wirepair<br>at start of freezing</td><td>Ω</td></tr><tr><td>α<sub>CW</sub></td><td>intensity parameter, CWconv</td><td>g·m·kg<sup>-1</sup>·°C<sup>-1</sup>·s<sup>-1</sup></td></tr><tr><td>α<sub>GN</sub></td><td>intensity parameter, GN13</td><td>kg·m<sup>-3</sup>·s<sup>-1</sup></td></tr><tr><td>α<sub>RJW</sub></td><td>intensity parameter, RJW14</td><td></td></tr><tr><td>α<sub>T</sub></td><td>intensity parameter, T13</td><td>m/s</td></tr><tr><td>β</td><td>saline density coefficient</td><td>kg/g</td></tr><tr><td>γ</td><td>electrical conductivity of brine</td><td>Ω<sup>-1</sup>·m<sup>-1</sup></td></tr><tr><td>γ<sub>0</sub></td><td>electrical conductivity of brine<br>at start of freezing</td><td>Ω<sup>-1</sup>·m<sup>-1</sup></td></tr><tr><td>κ</td><td>rhodamine adsorbed per<br>surface area</td><td>µg/cm<sup>2</sup></td></tr><tr><td>μ</td><td>brine dynamic viscosity</td><td>kg·m<sup>-1</sup>·s<sup>-1</sup></td></tr><tr><td>ν</td><td>brine kinematic viscosity</td><td>m<sup>2</sup>/s</td></tr></tbody></table>
<table><caption>**Table 1 Continued**</caption><thead><tr><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Π</td><td>local permeability</td><td>m<sup>2</sup></td></tr><tr><td>Π̄</td><td>effective permeability</td><td>m<sup>2</sup></td></tr><tr><td>ρ</td><td>brine density</td><td>kg/m<sup>3</sup></td></tr><tr><td>σ</td><td>reciprocal of the linear liquidus</td><td>°C·kg·g<sup>-1</sup></td></tr><tr><td>φ</td><td>liquid fraction</td><td></td></tr><tr><td>φ<sub>c</sub></td><td>critical liquid fraction</td><td></td></tr></tbody></table>

An alternative convective parameterization for gravity drainage was proposed by Griewank and Notz (2013). Their parameterization explicitly represents the upward and downward brine flow, assuming incompressible convective 2-D flow. The downwelling brine channel mass flux provides a direct connection between the sea ice interior and the ocean and is linearly dependent on the local Rayleigh number. Following Worster and Rees Jones (2015), we show in Appendix A that the original Griewank and Notz (2013) formulation can be recast into a vertical salt advection equation formally identical to equation (1), provided that the upwelling velocity in the kth layer is defined as

$$w_k = -\frac{\alpha_{\text{GN}}}{\rho} \sum_{i=1}^{k} (Ra_i - Ra_c) \Delta z_i. \quad (4)$$

We call this scheme “GN13.” Similar to RJW14, there are two tuning parameters: $\alpha_{GN}$, a convection strength tuning parameter, and $Ra_c$.

A third convective parameterization was proposed by Turner et al. (2013). In essence, it is intermediate between RJW14 and GN13 (Worster & Rees Jones, 2015). The potential upward flow rate is calculated in each of the model layers. It is then assumed that the actual flow rate is the largest among all model layers. The Turner et al. (2013) scheme also includes a term that relaxes the bulk salinity profile to some low-limiting value. The essential character of Turner et al. (2013) is given by

$$w_k = - \max_{z \le z_k} \alpha_T(z) (Ra(z) - Ra_c), \quad (5)$$

where $\alpha_T$ is a function capturing important properties of the mushy layer (Worster & Rees Jones, 2015). For each model layer, equation (5) sets $w$ as the largest value at or above that layer. We call this modified formulation of the Turner et al. (2013) parameterization “T13.”

## 2.2. A Convective Reformulation of Cox and Weeks (1988)

In a seminal study, Cox and Weeks (1988) proposed a gravity drainage parameterization derived from the measurements of Cox and Weeks (1975). This parameterization has since been widely used (Butler et al., 2016, 2017; Eicken, 1992; Saenz & Arrigo, 2012; Vancoppenolle et al., 2006, 2007). Their parameterization desalinates sea ice via two mechanisms. First, any newly formed sea ice is instantly desalinated by partitioning the salt between the new sea ice and ocean. Second, the sea ice desalinates by gravity drainage according to

$$\frac{\partial S}{\partial t} = -\alpha_{\text{CW}} \cdot \max \left( \frac{\phi}{\phi_c} - 1, 0 \right) \cdot \frac{\partial T}{\partial z}. \quad (6)$$

Equation (6) empirically relates the vertical sea ice temperature gradient, $\frac{\partial T}{\partial z}$, to the change in bulk salinity with time, $\frac{\partial S}{\partial t}$, using a desalination intensity parameter, $\alpha_{CW}$. Gravity drainage is moderated by the liquid fraction, $\phi$, and is switched off if $\phi$ falls below a critical liquid fraction threshold, $\phi_c$.

The Cox and Weeks (1988) parameterization is not suitable for this study for two reasons. First, equation (6) is incompatible with general solute transport because its empirical reliance on the temperature gradient may not be valid for solutes with a vertical brine concentration profile different to that of salinity. Second, segregation of salt at the advancing sea ice/ocean interface is incompatible with observations (Notz & Worster, 2008), if sea ice is defined to include the skeleton layer, as we do here.

We therefore reformulate the original Cox and Weeks (1988) parameterization as a convective scheme that is compatible with general solute transport. To do so, we first reject initial segregation, ensuring desalination
is solely by gravity drainage. We also recast equation (6) such that it can be solved alongside RJW14, GN13, and T13 using equation (1) by setting

$$w_{CW} = -\sigma \cdot \alpha_{CW} \cdot \max\left(\frac{\phi}{\phi_c} - 1, 0\right), \qquad (7)$$

where $\sigma = -0.0589 \, ^{\circ}\text{C}\cdot\text{kg}\cdot\text{g}^{-1}$ (Notz, 2005) is the reciprocal of the linear liquidus that converts the temperature gradient to a brine salinity such that $S_{br} = T/\sigma$. We keep $\phi_c = 0.05$ given the robustness of this sea ice permeability threshold in the literature (Golden et al., 1998; Pringle et al., 2009) and take $\alpha_{CW}$ to be a free tuning parameter. We call this convective reformulation of the Cox and Weeks (1988) parameterization “CWconv.”

## 2.3. Diffusive Parameterizations

Another stream of authors (Jeffery et al., 2011; Vancoppenolle et al., 2010) proposed representing gravity drainage with a diffusion equation (Figure 1)

$$\frac{\partial S}{\partial t} = \frac{\partial}{\partial z} \left( D_e \frac{\partial S_{br}}{\partial z} \right), \qquad (8)$$

where $D_e$ is an effective diffusivity. The physical mechanism behind diffusive schemes is not the molecular diffusion of salt which is much too weak to cause observed desalination rates (Notz & Worster, 2009; Untersteiner, 1968; Vancoppenolle et al., 2010). Rather, diffusive schemes are analogous to turbulent mixing in some ocean models. Convection mixes salt and tracers in the same way as molecular diffusion, but with stronger intensity. Intensified mixing is achieved by setting $D_e$ to a much larger value than molecular diffusivity.

Vancoppenolle et al. (2010) proposed a condition, from here called “VC10,” whereby $D_e$ switches from molecular to turbulent diffusion across a critical Rayleigh number threshold,

$$D_e(Ra) = \phi \cdot \{[1 - \chi(Ra)]D_m + \chi(Ra)D_{tur}\}, \qquad (9)$$

where a tunable reference turbulent diffusivity, $D_{tur}$, is set to around 3 orders of magnitude larger than the molecular diffusivity, $D_m$. The function $\chi(Ra)$ is based on a hyperbolic tangent function and acts as a smooth switch between $D_m$ and $D_{tur}$ across some critical Rayleigh number. When the liquid fraction decreases below a permeability threshold of $\phi_c = 0.05$ the sea ice is presumed to be impermeable and $D_e$ is set to 0.

One drawback of the VC10 formulation is that the diffusivity (and therefore the salt flux) would go to zero at the sea ice/ocean interface because $Ra$ vanishes for $z = h$. This effect is not present in results produced using VC10 (Kotovitch et al., 2016; Moreau et al., 2015; Vancoppenolle et al., 2010) because they computed the bottommost diffusivity based on the midpoint $Ra$. Hence, the VC10 numerical scheme does not exactly solve a diffusion equation and has an undesired dependence on vertical resolution.

Jeffery et al. (2011) proposed two alternative formulations for diffusivity (Jeffery & Hunke, 2014; Jeffery et al., 2011). The first, which they refer to as “enhanced molecular diffusivity,” sets

$$D_e = \phi D_{tur}. \qquad (10)$$

We call this parameterization “J11_EMD.” The other, which they call “mixing length diffusivity,” is to set the diffusivity as proportional to some mixing length

$$D_e = \begin{cases} \frac{g\Pi}{\mu} \Delta\rho l_m & \text{if } \Delta\rho \text{ is unstable} \\ 0 & \text{otherwise} \end{cases} \qquad (11)$$

where $\Pi$ is the local permeability, $\Delta\rho$ is the vertical brine density gradient (calculated using Jeffery et al., 2011, equation (16a)) and $l_m$ is a tunable mixing length parameter. We call this parameterization “J11_MLD.”

# 3. Methodology

## 3.1. Experiment Methodology

### 3.1.1. Laboratory Experiment Design

We performed a laboratory sea ice growth experiment to collect measurements against which to evaluate the parameterizations using the Roland von Glasow air-sea-ice chamber (RvG-ASIC). The RvG-ASIC is a
highly instrumented, open-topped glass tank (2.4 × 1.4 m² footprint, 1.2 m depth), contained in a cold room. The sides of the tank are insulated and can be gently heated from the outside. Using an appropriate level of heating, we are able to maintain free floating sea ice without creating visible open water at the tank sides. The tank is permanently open to a smaller side tank via four 10 cm holes. This side tank is never allowed to freeze over totally, and so provides a path for sample and spike lines into the ocean. Instrumentation to be frozen into the sea ice is mounted near the center of the tank on a vertical pole that is free to rise with the sea ice as it forms, and its freeboard increases (see also Thomas, 2018, Chapter 4).

In this experiment we ran separate sample lines into the tank for spiking rhodamine, and for extracting water for analysis. We mounted chains of digital thermometers at 1 cm vertical resolution to measure the in situ sea ice and ocean temperature (± 0.1 °C).

We made our artificial ocean using NaCl (AkzoNobel Sanal-P) and ultrapure water. The initial ocean salinity was (33.9 ± 1.1) g/kg (2σ analytical error). To dissolve the salt, and to maintain a well-mixed ocean, two pumps were mounted in opposite corners of the tank and ran at a low setting throughout the experiment. The pumps faced each other to minimize currents while generating turbulence to mix the water (Loose et al., 2011).

The ocean was cooled at an atmospheric temperature of -20 °C and the point of first sea ice formation was logged using a camera. During the experiment we changed the atmospheric temperature to -10 °C for 4 days before reverting back to -20 °C. In doing so we provided the sea ice with a reasonably complex forcing to give the models a difficult test. After around 10 cm of sea ice growth we injected 0.283 g of rhodamine 6G into the underlying water. At the end of the experiment the ocean rhodamine concentration was (89 ± 18) µg/kg (2σ), which is consistent with the mass of the system and the mass of the rhodamine spike. After around 17 cm of sea ice growth, 8.9 days since the first ice formation, we ended the experiment.

### 3.1.2. Slab Salinity and Rhodamine Profiles

Just after the experiment we collected two slabs of sea ice following protocols to minimize brine loss upon extraction (Cottier et al., 1999). We extracted the slabs in boxes while they were still floating on ocean, drained the ocean until the slabs were just touching the base of the boxes, and then shock froze the slabs with the ocean in a -40 °C freezer. Our only modification to the method of Cottier et al. (1999) was to place plastic sheets between the sea ice and ocean before shock freezing to prevent infiltration of salt and rhodamine from the freezing ocean into the sea ice sample. After several days of conditioning at -40 °C all of the liquid in the slabs was immobile and we removed the slabs from the freezer, taking them to a -25 °C cold room for processing. Using a bandsaw with a precleaned blade, we trimmed and discarded the now frozen ocean from the sides of the sample. We also removed at least 3 cm from the slabs at each edge to prevent any ocean contamination of the sample. The approximate dimensions of the trimmed slabs were ≈280 cm², 18 cm depth (slab 1); ≈500 cm², 17 cm depth (slab 2). Finally, we cut the trimmed slabs into 1 cm vertical layers.

Each layer was then melted and analyzed for salinity and rhodamine concentrations. The salinity was measured using a conductivity probe (CDC401 conductivity probe) against a set of volumetrically made standards. Rhodamine was measured by fluorometry (PerkinElmer LS45, 500 nm excitation, 525 nm emission) against a set of volumetrically made standards. The calculated density of the samples was used to convert the measured per liter concentrations to per kilogram. Analytical errors were calculated using standard analytical procedures (Taylor, 1997).

We rejected the bottommost layer from each slab from our analysis and our figures because there was an obvious artifact. The measured salinity in this layer was elevated over ocean salinity by 20% to 50% which we believe is due to contamination by the rejection of salt from the freezing ocean that surrounded the sample. Shallower layers were protected from this contamination in the -40 °C freezer by the plastic dividers.

### 3.1.3. Wireharp Salinity Profiles

We deployed wireharps to produce a 1 cm resolution vertical profile of liquid fraction (Notz & Worster, 2008; Notz et al., 2005). The liquid fraction was calculated using

$$ \phi = \frac{\gamma_0 Z_0}{\gamma(T, S_{br})Z} \qquad (12) $$

The $Z_0/Z$ term in equation (12) is the ratio of the resistance as the sea ice front just passes a given measuring depth to the resistance of that wirepair at some later point in time. We performed a baseline correction
and smoothed the measured resistance values over a 19 min window. On the advice of Dirk Notz we take $\frac{\gamma_0}{\gamma(T,S_{br})}$ to be equal to one. This term corrects for the increase in brine conductivity, $\gamma$, as the salinity increases relative to the conductivity at the start of freezing, $\gamma_0$. We cannot perform the necessary calibration for this term with our current electronics setup, but its effect on the resulting bulk salinity is expected to be small. $S$ is then recovered for each measurement depth using $S = \phi S_{br}(T)$. To calculate brine salinity, we assume thermodynamic equilibrium and use the liquidus curve derived by Rees Jones and Worster (2014) from the data of Weast (1971) for NaCl,

$$S_{br}/(\text{g kg}^{-1}) = -17.6 \left(\frac{T}{^{\circ}\text{C}}\right) - 0.389 \left(\frac{T}{^{\circ}\text{C}}\right)^2 - 0.00362 \left(\frac{T}{^{\circ}\text{C}}\right)^3. \quad (13)$$

The uncertainty on $S$ is calculated by Gaussian propagation of the uncertainty in the measured $T$ ($\pm 0.1^{\circ}\text{C}$) and estimated $Z_0$ ($\pm 0.5^{\circ}\text{C}$) (Zeigermann, 2018).

Rees Jones and Worster (2014) note that the wireharps systematically underestimate bulk salinity. We have found the same in previous experimental runs, particularly in cold, relatively impermeable sea ice. We do not understand this bias well but believe it relates to the increasing frequency, as liquid fraction decreases, of isolated brine inclusions which are less available to conduct current between the pairs of wires. We therefore restrict our use of the wireharps to the first 3.4 days of our experiment, where the liquid fraction as measured by the wireharps never drops below 0.9 and we presume there are few isolated brine inclusions. Restricting ourselves to this period also minimizes the relative error caused by neglecting the brine conductivity term, which will increase as brine salinity increases later in the experiment.

## 3.2. Model Methodology

### 3.2.1. Model Design

The model we use is based on two depth-dependent equations, one for salt (equation (1) or (8)), and the other for a second dissolved tracer, here rhodamine. The second dissolved tracer equation is formally similar to the salt equation, except that bulk salinity is replaced by bulk rhodamine concentration, $C$, and brine salinity is replaced by brine rhodamine concentration, $C_{br}$.

The model is forced by measured sea ice thickness and vertically resolved temperature profiles, linearly remapped onto a uniform grid with a fixed number of layers. We use measurements instead of calculations for $T$ and $h$ to minimize errors in our representation of the thermodynamics. To obtain the surface temperature $T_s$, we linearly extrapolated the vertical temperature profile to 0 depth. Finally, we calculate $h$ by linearly extrapolating the vertical temperature profile to the measured ocean temperature and smoothing the resulting thickness over a 10 min moving average. Our model is initialized with $h = 1$ cm. Each model layer initially has a bulk salinity equal to that of the underlying water at the start of the experiment, $S_w = 33.9$ g/kg.

$S_{br}$ is calculated using equation (13) with the exception of CWconv for which we use an inversion of the linear liquidus. This discrepancy is small at the encountered temperatures (Notz, 2005; Vancoppenolle et al., 2019). For liquid fraction, we neglect precipitated minerals, which is an excellent approximation for the encountered temperatures (Vancoppenolle et al., 2019), such that

$$\phi = S/S_{br}. \quad (14)$$

To get brine rhodamine concentration we assume the rhodamine is completely dissolved in brine (Appendix B: justifies this assumption) so that

$$C_{br} = \frac{C}{\phi}. \quad (15)$$

Finally, permeability is taken as a function of brine fraction (Freitag, 1999), with the effective permeability given by the harmonic mean approach (Griewank & Notz, 2015; Rees Jones & Worster, 2014).

### 3.2.2. Implementation of Parameterizations

We simulated salt and rhodamine vertical profiles using seven parameterizations, spanning the range of approaches in the literature (Table 2). We used four convective parameterizations in the framework of
<table><caption>**Table 2**

*Overview of the Gravity Drainage Parameterizations Evaluated in This Study*</caption><thead><tr><th>Abbreviation</th><th>Type</th><th>Tuning parameter</th><th>Tuned value</th><th>Reference value</th></tr></thead><tbody><tr><td>CWconv<sup>a</sup></td><td>Convective</td><td>α<sub>CW</sub> [g·m·(kg·°C·s)<sup>-1</sup>]</td><td>7.2 × 10<sup>-7</sup></td><td>n.a.</td></tr><tr><td>GN13<sup>b</sup></td><td>Convective</td><td>Ra<sub>c</sub>, α<sub>GN</sub> [kg·(m<sup>3</sup>·s)<sup>-1</sup>]</td><td>2.4, 6.7 × 10<sup>-3</sup></td><td>4.89, 5.84 × 10<sup>-4</sup></td></tr><tr><td>T13<sup>c</sup></td><td>Convective</td><td>Ra<sub>c</sub>, α<sub>T</sub> [m/s]</td><td>3.1, 6.4 × 10<sup>-7</sup></td><td>n.a.</td></tr><tr><td>RJW14<sup>d</sup></td><td>Convective</td><td>Ra<sub>c</sub>, α<sub>RJW</sub></td><td>2.9, 0.13</td><td>40, 0.03</td></tr><tr><td>VC10<sup>e</sup></td><td>Diffusive</td><td>Ra<sub>c</sub>, D<sub>tur</sub> [m<sup>2</sup>/s]</td><td>4.9, 5.5 × 10<sup>-8</sup></td><td>5, 10<sup>-6</sup></td></tr><tr><td>J11_EMD<sup>f</sup></td><td>Diffusive</td><td>D<sub>tur</sub> [m<sup>2</sup>/s]</td><td>2.8 × 10<sup>-8</sup></td><td>1.7 × 10<sup>-8</sup></td></tr><tr><td>J11_MLD<sup>g</sup></td><td>Diffusive</td><td>l<sub>m</sub> [m]</td><td>1.1 × 10<sup>-4</sup></td><td>7 × 10<sup>-5</sup></td></tr></tbody></table>

equation (1). Of these, RJW14, GN13, and T13 use the Rayleigh number to diagnose convection whereas CWconv uses $\phi$. We use CWconv rather than Cox and Weeks (1988) to enable us to simulate rhodamine and to remain consistent with the other parameterizations that reject initial segregation. For T13 we remove the slow relaxation term, which is incompatible with the transport of solutes other than salt, and take $\alpha_T$ to be a desalination strength tuning constant. The required boundary condition at the sea ice/ocean interface ($z = h$) is

$$S_{br} = S_{w}, \quad C_{br} = C_{w}. \qquad (16)$$

At $z = 0$ we take $S_{br} = S_{br}(T_s)$ (equation (13)) and $C_{br} = C_{br,s}$. Surface brine concentration of rhodamine is obtained from $C_{br,s} = C/\phi$, where $C$ and $\phi$ are obtained from linear extrapolation. For all four parameterizations, we used the same temporally explicit (first order in time), spatially centered (second order in space), conservative scheme, which was stable for the time step and grid layering used (25 layers). Stability was examined by looking at the CFL criterion.

We also include three diffusive parameterizations (Table 2). The first is VC10, retained in spite of its drawbacks because it was used in several biogeochemical studies (Kotovitch et al., 2016; Moreau et al., 2015; Vancoppenolle et al., 2010). We also include J11_EMD and J11_MLD (Jeffery et al., 2011). The boundary conditions we used for diffusive schemes are

$$\begin{cases} F_s = 0, & F_c = 0 & \text{for } z = 0, \\ S_{br} = S_w, & C_{br} = C_w & \text{for } z = h, \end{cases} \qquad (17)$$

where $F_s = -D_e \cdot \partial S_{br}/\partial z$ and $F_c = -D_e \cdot \partial C_{br}/\partial z$ are the diffusive salt and rhodamine fluxes, respectively. For diffusive schemes, we use a temporally implicit (first order in time), spatially centered (second order in space), conservative, unconditionally stable scheme (Vancoppenolle et al., 2010). All other details are provided with the model code (see Acknowledgments section for repository).

### 3.2.3. Tuning Parameters

To treat all parameterizations fairly we tune each to the same target. We choose the mean of the measured $S$ of slabs 1 and 2 as our target (Figure 2, bottom panel). The modeled $S$ is linearly mapped onto the mean slab $S$ using normalized depth, $z/h$, as the vertical coordinate. We use normalized depth because of the 1 cm difference in slab thickness.

The performance of the tuning parameters is quantified using the mean absolute deviation ($d_{MA}$) of the modeled ($S^{\text{mod}}$) and measured ($S^{\text{obs}}$) bulk salinity for $N$ measured layers.

$$d_{MA} = \frac{\sum_k |S_k^{\text{mod}} - S_k^{\text{obs}}|}{N} \qquad (18)$$

The best tuning parameter or parameter pair (Table 2) was identified as that which produced the minimum $d_{MA}$.

<footnote><sup>a</sup>Cox and Weeks (1988), modified. <sup>b</sup>Griewank and Notz (2013) with reference parameter values from Griewank and Notz (2015). <sup>c</sup>Turner et al. (2013), modified. <sup>d</sup>Rees Jones and Worster (2014). <sup>e</sup>Vancoppenolle et al. (2010).</footnote>

<footnote><sup>f</sup>Enhanced-molecular diffusivity (EMD) parameterization from Jeffery et al. (2011). <sup>g</sup>Mixing-length diffusivity (MLD) parameterization from Jeffery et al. (2011), with reference parameter values from Jeffery and Hunke (2014).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Observed (symbols and black line) versus simulated (convective as red, diffusive as blue) vertical bulk salinity profiles. (top) Wireharp measurements early in the experiment. (bottom) Slab-based values at the end of the experiment. Vertical bars show the depth of sea ice that the measured layers integrate. Horizontal bars represent $2\sigma$ analytical error.</figcaption></figure>

There are some large differences between our calculated tuning parameters and those presented in the original papers (Table 2). Our lower value for $D_{\text{tur}}$ for VC10 is likely due to different time scales for the tuning target (several weeks for (Vancoppenolle et al., 2010), compared to several days for this study). Our $\alpha_{\text{GN}}$ is larger than those presented in Griewank and Notz (2013) and Griewank and Notz (2015). However, inspecting our parameter space we find low sensitivity to this parameter, and similar results could have been achieved using the original parameter (see supplementary information). We have a much lower $Ra_c$ when compared to those presented in Rees Jones and Worster (2014). Using an $Ra_c$ of 20 or 40 our simulated desalination was too weak. Rees Jones and Worster (2014) constrained their tuning parameters partially by enforcing a delay in the onset of gravity drainage, as observed in sea ice grown from a cold plate (Wettlaufer et al., 1997) but not in the field (Notz & Worster, 2008). Rees Jones and Worster (2014) note that without this constraint, as in this study, much lower values of $Ra_c$ are potentially reasonable.

**4. Results**

The wireharp and slab bulk salinity profiles are “C” shaped, with increased salinity near the sea ice base and the sea ice atmosphere interface (Figure 2). The salinity of the interior sea ice decreases through time. There is excellent agreement between the two slab profiles for salinity except for the shallowest measurement
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Observed (symbols and black line) versus simulated (convective as red, diffusive as blue) vertical bulk rhodamine concentration profiles, retrieved from slab samples taken at the end of the experiment. Vertical bars show the depth of sea ice that the measured layers integrate. Horizontal bars represent $2\sigma$ analytical error.</figcaption></figure>

depth. The slab rhodamine profiles also show increased rhodamine concentrations toward the sea ice base (Figure 3). Above this, between 0.9 and 0.5 normalized depth, is a region of near-homogeneous rhodamine concentration. The rhodamine concentration then begins to fall in the sea ice above around 0.5 normalized depth. Near-surface rhodamine concentrations are very low. As with salinity, there is excellent agreement between the rhodamine profiles in the two slabs.

Qualitatively, all of the parameterizations capture the increase in salinity near the sea ice base and the decrease in interior salinity through time. All of the schemes also capture the increase in rhodamine at around 0.9 normalized depth and the decrease in rhodamine toward the upper interface. There are, however, significant qualitative differences between the parameterizations.

Early on in the experiment there is a large spread in the predicted bulk salinities between the parameterizations (Figure 2). VC10 desalinates too much in the upper portion of the sea ice and overestimates salinity near the base. CWconv underestimates desalination for the entire 0.5 day profile. RJW14, GN13, and T13 overestimate desalination in the interior during the first 2 days of the experiment. The parameterizations also differ near the sea ice surface. RJW14, GN13, and J11_MLD consistently predict a “C”-shaped profile while VC10 and J11_EMD predict minimum salinities near the surface. CWconv predicts a “C” shape early on but shows no large increase in near-surface salinity at the end of the experiment, while T13 predicts only a small increase salinity near the top of the sea ice. The increase in near-surface salinity predicted by RJW14 and GN13 is confined to the top centimeter of the sea ice, which is consistent with the observations. In contrast, J11_MLD predicts an increase in near-surface salinity beginning toward the middle of the sea ice. Away from the upper interface, the four convective parameterizations perform similarly for the 3 day profile, and the differences between these schemes are larger earlier in the experiment. RJW14 and GN13 are qualitatively similar throughout.

With regard to rhodamine (Figure 3), the region of near homogenous concentration is not captured by the diffusive schemes, or by CWconv. These schemes tend to underestimate rhodamine concentrations between 0.3 and 0.7 normalized depth. The Rayleigh number based convective schemes, RJW14, GN13, and T13, capture the shape of the measured profile, though the depth at which rhodamine begins to fall at around 0.5 normalized depth is underestimated. RJW14, GN13, and T13 are remarkably consistent with each other when predicting rhodamine concentration.

We now turn to a quantitative evaluation of the parameterizations. To do so we calculate the residuals between observed and modeled concentrations after linearly remapping the modeled concentrations onto the observation depths. We then calculate the mean bias, $b_M$, and mean absolute deviation, $d_{MA}$, and use these as performance metrics (Table 3). For the wireharps we calculate $b_M$ and $d_{MA}$ for each measured profile during the period of retained data (3.4 days) and use the mean of $b_M$ and $d_{MA}$ over this period.

RJW14 or GN13 perform best for five of the six observational metrics, displaying the lowest $b_M$ and $d_{MA}$ for wireharp salinity and slab rhodamine (Table 3). T13 produces the joint best $d_{MA}$ for wireharp salinity and CWconv produces the worst $d_{MA}$ for slab salinity. J11_MLD achieves the lowest $b_M$ for slab salinity, though this reflects compensating regions of large over and underestimation (Figure 2). VC10 performs worst against the wireharps and J11_MLD performs worst against slab rhodamine. The performance of J11_EMD is always intermediate.

## 5. Discussion

We tuned all parameterizations against the same target (slab salinity), then evaluated them against this target and two different observational data sets (wireharp salinity early in the experiment and slab rhodamine). The performance of the parameterizations against slab salinity tells us how well these schemes can perform
<table><caption>**Table 3**

Evaluation Statistics of Simulated Bulk Salinity (S) and Bulk Rhodamine Concentration (C) by Comparison With Slab and Wireharp Observational Values</caption><thead><tr><th rowspan="2"></th><th colspan="2">Wireharp S</th><th colspan="2">Slab S</th><th colspan="2">Slab C</th></tr><tr><th></th><th>/ (g/kg)</th><th></th><th>/ (g/kg)</th><th></th><th>/ (µg/kg)</th></tr><tr><th></th><th>b<sub>M</sub></th><th>d<sub>MA</sub></th><th>b<sub>M</sub></th><th>d<sub>MA</sub></th><th>b<sub>M</sub></th><th>d<sub>MA</sub></th></tr></thead><tbody><tr><td>CWconv</td><td>-1.3</td><td>3.6</td><td>-1.5</td><td>3.1</td><td>0.6</td><td>3.9</td></tr><tr><td>GN13</td><td>0.1</td><td>2.6</td><td>-0.6</td><td>1.7</td><td>0.4</td><td>3.1</td></tr><tr><td>RJW14</td><td>-0.2</td><td>2.4</td><td>-0.6</td><td>1.4</td><td>-0.2</td><td>2.1</td></tr><tr><td>T13</td><td>-1.2</td><td>2.4</td><td>-0.8</td><td>1.9</td><td>0.6</td><td>2.3</td></tr><tr><td>VC10</td><td>-2.1</td><td>6.8</td><td>-2.6</td><td>2.7</td><td>-2.6</td><td>2.7</td></tr><tr><td>J11_EMD</td><td>-1.2</td><td>3.2</td><td>-0.9</td><td>2.8</td><td>-1.4</td><td>3.5</td></tr><tr><td>J11_MLD</td><td>2.1</td><td>3.5</td><td>0.2</td><td>2.9</td><td>-4.3</td><td>4.3</td></tr></tbody></table>

Wireharp data cover the first 3.4 days of the experiment. The evaluation statistics are the mean bias, $b_M$, and the mean absolute deviation, $d_{MA}$, between the measurements and model.

when predicting bulk sea ice salinity. Their performance against the wireharps gives us insights into the robustness of their performance as salt dynamics parameterizations. Their performance against rhodamine, with a vertical brine concentration profile deliberately decoupled from the salinity, tells us how well these salt dynamics parameterizations perform as brine dynamics parameterizations, and hence how useful they are when predicting the dynamics of any solute in growing sea ice.

Our experimental system is a reasonable approximation of natural sea ice with respect to gravity drainage. The initial ocean salinity and variation in experimental air temperatures are within the range of those naturally found during the sea ice growth season in the Arctic and Antarctic. We have also maintained free floating sea ice and so have simulated a natural freeboard. Doing so prevented flooding of the sea ice surface and, given that our freeboard was < 2 cm, we expect the freeboard to have had minimal impact on our results. In nature, growing sea ice always has a region adjacent to the ocean that exchanges solutes. In our experiments this region is not fundamentally different to similar regions in thicker, growing sea ice. The conclusions drawn from this study are therefore relevant for natural sea ice undergoing or having recently undergone gravity drainage. Our experiments are, however, shorter than the sea ice growth season and have simple forcing relative to the environment. Longer experiments with more realistic forcing could highlight additional differences between the schemes and may require further parameter adjustments.

Our modeling assumes that the rhodamine is perfectly dissolved. Part of the reason we chose rhodamine as a tracer is that the behavior of rhodamine in brine and at icy surfaces has been experimentally determined (Krembs et al., 2000). In the experiments of Krembs et al. (2000) rhodamine was soluble in cold and highly saline brines at the maximum concentrations in our experimental system. Rhodamine was, however, observed to adsorb to icy surfaces. To assess the impact of rhodamine adsorption in our study, we added an adsorption scheme to the model based on the experiments of Krembs et al. (2000) (Appendix B:). The effect of this adsorption does not affect the shape of the modeled bulk rhodamine profile and changes the bulk concentrations by less than 10%. Rhodamine is therefore appropriate as a dissolved tracer in this study.

Our formulations of RJW14, VC10, J11_EMD, and J11_MLD only differ from their original formulations in our choice of tuning parameters. Our reformulation of GN13 is practically identical to the version presented in Griewank and Notz (2013) and Griewank and Notz (2015) (Figure A1), with revised tuning parameters. There are, however, important differences between CWconv and T13 compared to the parameterization originally presented (Cox & Weeks, 1988; Turner et al., 2013). For CWconv the most significant of these differences is our rejection of initial segregation of salt at the advancing interface which was a key mechanism in Cox and Weeks (1988). The effects of this change can be seen in Figure A1. We therefore do not claim to evaluate the Cox and Weeks (1988) parameterization or the usefulness of the tools they developed. Rather, we use CWconv to investigate the effect of the increased complexity in RJW14, GN13, and T13. CWconv is an extremely simple gravity drainage parameterization, with the intensity of desalination mediated by the
vertical temperature (and, by extension, brine density) gradient and the local liquid fraction. RJW14, GN13, and T13 claim improved physical realism by replacing the liquid fraction with a more complex variable, the Rayleigh number. The difference between the performance of CWconv and the Rayleigh number based convective parameterizations is an estimate of the payoff from this increased complexity. T13 can be viewed as the simplest version of the Turner et al. (2013) parameterization, with complex terms wrapped up in a single tuning parameter, $\alpha_T$, and no slow desalination mode. The performance of T13 therefore (1) represents a lower bound on the potential performance of the Turner et al. (2013) parameterization and (2) explores the potential for using Turner et al. (2013) to transport solutes other than salt.

All of the parameterizations used in this study do a reasonable job in accurately predicting the final slab salinity, with absolute mean biases ranging from 0.2 (J11_MLD) to 2.6 (VC10) g/kg. All of the parameterizations can be useful tools when predicting the bulk salinity of young sea ice. When compared to salinity during the first 3.4 experimental days the parameterizations produce similar accuracy, with absolute mean biases ranging from 0.1 (GN13) to 2.1 (VC10) g/kg. The performance of these parameterizations therefore appears to be robust over the experimental period. However, inspecting the bulk salinity profiles during this period (Figure 2) we find that VC10 overestimates desalination in the upper portion of the sea ice and produces salinities in excess of the ocean near the sea ice base. The performance of VC10 when predicting salinity is less robust over the duration of our experiment than the other parameterizations. We also see differences between RJW14 and GN13 early in the experiment, particularly for the 0.5 day profile. While these parameterizations converge to very similar solutions they take slightly different routes. Supplementary information Figure S1 provides a timeseries of integrated sea-ice bulk salinity. T13 performs similarly to RJW14 and GN13 when predicting slab salinity and is around 1 g/kg less accurate when predicting wireharp salinity. RJW14 and GN13 are more accurate than for the simpler convective parameterization, CWconv, by around 1 g/kg for both observational metrics for salinity. RJW14 and GN13 perform similarly for all salinity metrics.

The performance of all the parameterizations against rhodamine concentrations is qualitatively reasonable. J11_MLD ($b_M$ = -4.3 µg/kg) underestimates rhodamine concentrations more than the other parameterizations (range of $b_M$ 0.6 to -2.6 µg/kg), which probably reflects the underestimation of bulk salinity in the lower half of the sea ice. The diffusive schemes underestimate rhodamine concentrations. The convective schemes are more accurate than the diffusive schemes (absolute mean biases of 0.2 to 0.6 and 1.4 to 4.3 µg/kg, respectively). As with salinity, RJW14 and GN13 (absolute mean biases of 0.2 and 0.4 µg/kg, respectively), are more accurate than the simpler CWconv (absolute mean bias 0.6 µg/kg). T13 performs similarly to RJW14 and GN13, highlighting the potential of T13 for transporting solutes other than salt. The outstanding accuracy of RJW14, GN13, and T13 when predicting rhodamine concentration is consistent with the analytical work of Rees Jones and Worster (2013b), and allows us to confidently use RJW14, GN13, and T13 as brine dynamics parameterizations, useful for predicting the dynamics of any solute.

RJW14 performed best overall, with the best or joint best performance for five of the six metrics, though the performance of RJW14 was similar to GN13 and T13. We use $d_{MA}$ (Table 3) as an estimate of the precision of RJW14. The mean bulk salinity measured by the wireharps over the first 3.4 days is 13 g/kg and the mean of the slab bulk salinity is 8.6 g/kg. RJW14 produces $d_{MA}$ of 2.4 g/kg (wireharp salinity) and 1.4 g/kg (slab salinity), resulting in precisions of 18% and 16%, respectively. The mean slab bulk rhodamine concentration is 12.4 µg/kg and the $d_{MA}$ for RJW14 is 2.1 µg/kg, producing a precision of 17%. We estimate the precision of RJW14 when predicting the bulk concentration profile of any solute in young or recently convecting sea ice is less than 20%.

## 6. Conclusions

We have evaluated seven gravity drainage parameterizations against tracer measurements in growing sea ice. RJW14 (Rees Jones & Worster, 2014), GN13 (Griewank & Notz, 2013), and a simplified version of Turner et al. (2013) (T13), which are based explicitly on convection diagnosed using a Rayleigh number, performed best, generally outperforming diffusive approaches and a simpler convective parameterization. RJW14, GN13, and T13 did an excellent job against salinity measurements, confirming that these parameterizations are powerful tools for predicting salt dynamics in growing sea ice. RJW14, GN13, and T13 performed similarly well for a dissolved tracer decoupled from the salinity, showing that they can be thought of as brine dynamics parameterizations and are powerful tools for predicting the dynamics of any solute in growing sea ice.
## Appendix: A: Equivalence of Convective Parameterizations

Worster and Rees Jones (2015) argue that RJW14 and GN13 are essentially similar in terms of underlying physics. Here we show that RJW14 and GN13, and indeed T13, can be expressed in the same numerical framework.

Let us first show that GN13 can be recast into the framework of RJW14. The GN13 salinity scheme is based on the local conservation of brine mass in the kth layer between upwelling (in the mush, $m_{br,k\uparrow}$) and downwelling (through brine channels, $m_{br,k\downarrow}$), assuming incompressible flow

$$m_{br,k\uparrow} = m_{br,k-1\uparrow} + m_{br,k\downarrow} = \sum_{i=1}^{k} m_{br,i\downarrow}. \quad (A1)$$

In this framework, the conservation of salt reads

$$\frac{S_k^{n+1} - S_k^n}{\Delta t} = \frac{m_{br,k\uparrow}}{\rho \Delta t} \left( \frac{S_{k+1/2}^{br,n} - S_{k-1/2}^{br,n}}{\Delta z} \right). \quad (A2)$$

Equation (A2) is equivalent to the RJW14 salt advection equation (equation (1)) provided

$$w_k = m_{br,k\uparrow} / (\rho \Delta t). \quad (A3)$$

GN13 impose that the brine downwelling mass is proportional to the local Rayleigh number

$$m_{br,k\downarrow} = \alpha_{GN}(Ra_k - Ra_c)\Delta z_k \Delta t. \quad (A4)$$

Inserting equation (A1) into equation (A3), and using A4 to reexpress $m_{br,i\downarrow}$, we get

$$w_k = -\frac{\alpha_{GN}}{\rho} \sum_{i=1}^{k} (Ra_i - Ra_c) \Delta z_i. \quad (A5)$$

<figure><img src="image_8.png" /><figcaption>**Figure A1.** Comparison of original formulations to convective reformulations for Cox and Weeks (1988) and CWconv (left), and Griewank and Notz (2013) expressed in its original form, and recast in the RJW14 framework, which we call GN13 (right).</figcaption></figure>
which expresses the RJW14 upwelling velocity as a function of the GN13 parameters. Our recast version of GN13 is virtually identical to the original (Figure A1).

The reverse operation, using the RJW14 upwelling velocity in the GN13 framework, is easier. One just needs to define

$$m_{br,k\uparrow} = \rho w_k \Delta t \qquad (A6)$$

$$m_{br,k\downarrow} = (w_k - w_{k+1})\rho\Delta t \qquad (A7)$$

Hence both schemes can be solved using the same framework. The main difference is the handling of the downwelling brine flux and vertical velocity.

Our reformulation of the parameterization proposed by Cox and Weeks (1988), CWconv, differs significantly from the original because we have omitted initial segregation of salt at the advancing sea ice interface (Figure A1).

## Appendix B:: Sensitivity to Rhodamine Adsorption

The adsorption of rhodamine 6G to icy surfaces was noted by Krembs et al. (2000). In our modeling, we assume that rhodamine is perfectly dissolved. We therefore added two adsorption schemes to our model to test the sensitivity of our results to this adsorption.

The first scheme uses equations from Krembs et al. (2000) to calculate the rhodamine adsorbed to the internal surfaces of the sea ice. This scheme uses empirical fits relating rhodamine adsorption and internal sea ice area (Krembs et al., 2000) to calculate the amount of rhodamine adsorbed. The internal sea ice area, A, is parameterized using

$$A/(\text{cm}^2/\text{kg}) = e^{0.63(T/^{\circ}\text{C})} \times 1.35 \times 10^5. \qquad (B1)$$

Equation (B1) was derived from data for laboratory grown columnar sea ice having undergone some desalination. The rhodamine 6G adsorbed per surface area, κ, is calculated using

$$\kappa/(\mu\text{g/cm}^2) = 11.413 \times (C_{\text{br}}/(\mu\text{g/kg})) - 1.7268 \times 10^{-4}, \qquad (B2)$$

<figure><img src="image_9.png" /><figcaption>**Figure B1.** Vertical profiles of rhodamine measured in the two slab samples and modeled concentrations of a rhodamine at the end of the experiment. In all cases the model is run using RJW14. The model output with no adsorption is shown alongside the Krembs and simple adsorption schemes. The vertical bars show the depth of sea ice that the measured layers integrate. The horizontal bars represent 2σ analytical error.</figcaption></figure>
where $C_{br}$ is the rhodamine concentration in the brine. The rhodamine adsorbed to the internal sea ice surfaces, $C_{ads}$, is then given by

$$C_{ads} = \kappa \cdot A. \qquad (B3)$$

We call this scheme “Krembs.”

We also added a second, simpler adsorption scheme to the model. In this scheme

$$C_{ads} = x_{ads}C_{br}, \qquad (B4)$$

where $x_{ads}$ is a constant between 0 and 1 that partitions the rhodamine between an adsorbed and aqueous phase. We call this second adsorption scheme “Simple.”

$C_{br}$ after adsorption is diagnosed by mass conservation. This reduced $C_{br}$ is then transported by gravity drainage while $C_{ads}$ remains stationary. The bulk concentration of rhodamine, C, is recovered using

$$C = (1 - \phi)(C_{ads} + C_{br}) \qquad (B5)$$

The impact of the two adsorption schemes is not dramatic (Figure B1). The shape of the tracer profile is preserved for both adsorption schemes. For Krembs, which we believe is the most realistic, the maximum difference from the base case is around 10%. Rhodamine 6G is a reasonable approximation of a perfectly dissolved solute for the purposes of this study.

## Acknowledgments

This article is dedicated to the memory of Professor Roland von Glasow (30.07.1971 – 06.09.2015), who was the inspiration and driving force behind the development of the air-sea-ice chamber that is described here. We would like to thank Stuart Rix, Dave Blomfield, Nick Griffin, Andy Macdonald, and Gareth Flowerdue for technical support with the RvG-ASIC. Thanks also to Dirk Notz and Leif Riemenschneider for many stimulating discussions and technical help with the wireharps, to David Rees Jones and Grae Worster for assistance in setting up their parameterization, and to two anonymous reviewers for their thorough evaluations of the manuscript. This work was carried out under the auspices of BEPSII (Biogeochemical Exchange Processes at the Sea-Ice Interfaces), a network supported by the Scientific Committee on Antarctic Research (SCAR), the Climate and Cryosphere (CliC) project of the World Climate Research Programme (WCRP) and by the Surface Ocean Lower Atmosphere Study (SOLAS). Dorothee Bakker was funded by the Natural Environment Research Council (NERC) PICCOLO (Processes Influencing Carbon Cycling: Observations of the Lower limb of the Antarctic Overturning) project (NE/P021395/1). This work was supported by project EISPAC (NERC Grant no. NE/R012857/1), part of the Changing Arctic Ocean programme, jointly funded by the UKRI Natural Environment Research Council (NERC) and the German Federal Ministry of Education and Research (BMBF). This work received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7-2007-2013, Grant Agreement 616938) and the Horizon 2020 research and innovation program through the EUROCHAMP-2020 Infrastructure Activity under Grant Agreement 730997. After R. v. G.'s death the School of Environmental Sciences at UEA kindly funded M. T.'s PhD stipend. All data, plot scripts, and model code are available at https://doi.org/10.5281/ zenodo.3505780 (submission version) and https://doi.org/10.5281/zenodo. 3626975 (final version).

## References

Assur, A. (1958). Composition of sea ice and its tensile strength. *Arctic Sea Ice*, 598, 106–138.

Butler, B. M., Papadimitriou, S., Day, S. J., & Kennedy, H. (2017). Gypsum and hydrohalite dynamics in sea ice brines. *Geochimica et Cosmochimica Acta*, 213, 17–34.

Butler, B. M., Papadimitriou, S., Santoro, A., & Kennedy, H. (2016). Mirabilite solubility in equilibrium sea ice brines. *Geochimica et Cosmochimica Acta*, 182, 40–54.

Cottier, F., Eicken, H., & Wadhams, P. (1999). Linkages between salinity and brine channel distribution in young sea ice. *Journal of Geophysical Research*, 104(C7), 15,859–15,871.

Cox, G. F., & Weeks, W. F. (1975). Brine drainage and initial salt entrapment in sodium chloride ice. (Tech. Rep.): DTIC Document.

Cox, G. F., & Weeks, W. F. (1988). Numerical simulations of the profile properties of undeformed first-year sea ice during the growth season. *Journal of Geophysical Research*, 93(C10), 12,449–12,460. https://doi.org/10.1029/JC093iC10p12449

Eicken, H. (1992). Salinity profiles of Antarctic sea ice: Field data and model results. *Journal of Geophysical Research*, 97(C10), 15,545–15,557.

Eide, L., & Martin, S. (1975). The formation of brine drainage features in young sea ice. *Journal of Glaciology*, 14(70), 137–154.

Feltham, D., Untersteiner, N., Wettlaufer, J., & Worster, M. (2006). Sea ice is a mushy layer. *Geophysical Research Letters*, 33, L14501. https:// doi.org/10.1029/2006GL026290

Freitag, J. (1999). Untersuchungen zur hydrologie des arktischen meereises: Konsequenzen für den kleinskaligen stofftransport (in German) the hydraulic properties of Arctic sea-ice: Implications for the small scale particle transport. *Reports on polar research*, 325.

Fripiat, F., Meiners, K. M., Vancoppenolle, M., Ackley, S. F., Arrigo, K. R., Carnat, G., et al. (2017). Macro-nutrient concentrations in Antarctic pack ice: Overall patterns and overlooked processes. *Elementa*, 5–13, 1–24.

Fritsen, C., Lytle, V., Ackley, S., & Sullivan, C. (1994). Autumn bloom of Antarctic pack-ice algae. *Science*, 266(5186), 782–784.

Golden, K., Ackley, S., & Lytle, V. (1998). The percolation phase transition in sea ice. *Science*, 282(5397), 2238–2241.

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, 3370–3386. https://doi.org/10.1002/jgrc.20247

Griewank, P. J., & Notz, D. (2015). A 1-D modelling study of Arctic sea-ice salinity. *The Cryosphere*, 9(1), 305–329.

Jeffery, N., & Hunke, E. (2014). Modeling the winter-spring transition of first-year ice in the Western Weddell Sea. *Journal of Geophysical Research: Oceans*, 119, 5891–5920. https://doi.org/10.1002/2013JC009634

Jeffery, N., Hunke, E., & Elliott, S. (2011). Modeling the transport of passive tracers in sea ice. *Journal of Geophysical Research*, 116, C07020. https://doi.org/10.1029/2010JC006527

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 Part II: Topical Studies in Oceanography*, 131, 75–83.

König, D., Miller, L. A., Simpson, K. G., & Vagle, S. (2018). Carbon dynamics during the formation of sea ice at different growth rates. *Frontiers in Earth Sciences*, 6(234). https://doi.org/10.3389/feart.2018.00234

Kotovitch, M., Moreau, S., Zhou, J., Vancoppenolle, M., Dieckmann, G. S., Evers, K.-U., et al. (2016). Air-ice carbon pathways inferred from a sea ice tank experiment. *Elementa Science of the Anthropocene*, 4, 000112.

Krembs, C., Gradinger, R., & Spindler, M. (2000). Implications of brine channel geometry and surface area for the interaction of sympagic organisms in Arctic sea ice. *Journal of Experimental Marine Biology and Ecology*, 243(1), 55–80.

Lake, R., & Lewis, E. (1970). Salt rejection by sea ice during growth. *Journal of Geophysical Research*, 75(3), 583–597.

Loose, B., Miller, L. A., Elliott, S., & Papakyriakou, T. (2011). Sea ice biogeochemistry and material transport across the frozen interface. *Oceanography*, 24(3), 202–218.

Middleton, C., Thomas, C., De Wit, A., & Tison, J.-L. (2016). Visualizing brine channel development and convective processes during artificial sea-ice growth using Schlieren optical methods. *Journal of Glaciology*, 62(231), 1–17.

Moreau, S., Vancoppenolle, M., Delille, B., Tison, J.-L., Zhou, J., Kotovitch, M., et al. (2015). Drivers of inorganic carbon dynamics in first-year sea ice: A model study. *Journal of Geophysical Research: Oceans*, 120, 471–495. https://doi.org/10.1002/2014JC010388
Moreau, S., Vancoppenolle, M., Zhou, J., Tison, J.-L., Delille, B., & Goosse, H. (2014). Modelling argon dynamics in first-year sea ice. *Ocean Modelling*, 73, 1–18. https://doi.org/10.1016/j.ocemod.2013.10.004

Niedrauer, T. M., & Martin, S. (1979). An experimental study of brine drainage and convection in young sea ice. *Journal of Geophysical Research*, 84(C3), 1176–1186.

Notz, D. (2005). Thermodynamic and fluid-dynamical processes in sea ice (PhD thesis).

Notz, D., Wettlaufer, J. S., & Worster, M. G. (2005). A non-destructive method for measuring the salinity and solid fraction of growing sea ice in situ. *Journal of Glaciology*, 51(172), 159–166.

Notz, D., & Worster, M. G. (2008). In situ measurements of the evolution of young sea ice. *Journal of Geophysical Research*, 113, C03001. https://doi.org/10.1029/2007JC004333

Notz, D., & Worster, M. G. (2009). Desalination processes of sea ice revisited. *Journal of Geophysical Research*, 114, C05006. https://doi.org/10.1029/2008JC004885

Pellichero, V., Sallée, J.-B., Chapman, C. C., & Downes, S. M. (2018). The Southern Ocean meridional overturning in the sea-ice sector is driven by freshwater fluxes. *Nature Communications*, 9, 1789.

Pringle, D., Miner, J., Eicken, H., & Golden, K. (2009). Pore space percolation in sea ice single crystals. *Journal of Geophysical Research*, 114, C12017. https://doi.org/10.1029/2008JC005145

Rees Jones, D. W., & Worster, M. G. (2013a). Fluxes through steady chimneys in a mushy layer during binary alloy solidification. *Journal of Fluid Mechanics*, 714, 127–151.

Rees Jones, D. W., & Worster, M. G. (2013b). A simple dynamical model for gravity drainage of brine from growing sea ice. *Geophysical Research Letters*, 40, 307–311. https://doi.org/10.1029/2012GL054301

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, 5599–5621. https://doi.org/10.1002/2013JC009296

Saenz, B. T., & Arrigo, K. R. (2012). Simulation of a sea ice ecosystem using a hybrid model for slush layer desalination. *Journal of Geophysical Research*, 117, C05007. https://doi.org/10.1029/2011JC007544

Taylor, J. (1997). *Introduction to error analysis, the study of uncertainties in physical measurements*. Sausalito, Calif: University Science Books.

Thomas, M. (2018). Brine and pressure dynamics in growing sea ice: Measurements and modelling in the Roland von Glasow air-sea-ice chamber (PhD thesis), University of East Anglia.

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*, 118(5), 2279–2294. https://doi.org/10.1002/jgrc.20171

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

Vancoppenolle, M., Bitz, C. M., & Fichefet, T. (2007). Summer landfast sea ice desalination at Point Barrow, Alaska: Modeling and observations. *Journal of Geophysical Research*, 112, C04022. https://doi.org/10.1029/2006JC003493

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

Vancoppenolle, M., Goosse, H., De Montety, A., Fichefet, T., Tremblay, B., & Tison, J.-L. (2010). Modeling brine and nutrient dynamics in Antarctic sea ice: The case of dissolved silica. *Journal of Geophysical Research*, 115, C02005. https://doi.org/10.1029/2009JC005369

Vancoppenolle, M., Madec, G., Thomas, M., & McDougall, T. J. (2019). Thermodynamics of sea ice phase composition revisited. *Journal of Geophysical Research: Oceans*, 124, 615–634. https://doi.org/10.1029/2018JC014611

Vancoppenolle, M., Meiners, K. M., Michel, C., Bopp, L., Brabant, F., Carnat, G., et al. (2013). Role of sea ice in global biogeochemical cycles: Emerging views and challenges. *Quaternary Science Reviews*, 79, 207–230.

Weast, R. C. (1971). *Handbook of chemistry and physics* (52nd ed.): Chemical Rubber Co.

Wells, A., Wettlaufer, J., & Orszag, S. (2011). Brine fluxes from growing sea ice. *Geophysical Research Letters*, 38, L04501. https://doi.org/10.1029/2010GL046288

Wettlaufer, J., Worster, M. G., & Huppert, H. E. (1997). The phase evolution of young sea ice. *Geophysical Research Letters*, 24(10), 1251–1254.

Worster, M. G. (1992). *Interactive dynamics of convection and solidification* *Edited by* Davis, S. H., Huppert, H. E., Müller, U., & Worster, M. G., pp. 113–138. Netherlands: Kluwer.

Worster, M. G., & Rees Jones, D. W. (2015). Sea-ice thermodynamics and brine drainage. *Philosophical Transactions Royal Society A*, 373(2045), 20140166.

Zeigermann, L. M. (2018). Desalination processes in thin sea ice during freezing and melting (Master's Thesis), University Hamburg, Hamburg, Germany.