Tracer Measurements in Growing Sea Ice Support Convective Gravity Drainage Parameterizations

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. 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.


Introduction
Vast areas of the world's oceans freeze each year. Some of the salt originally contained in the freezing seawater 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 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., 2014Moreau et al., , 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.

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 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. of them, thermodynamic equilibrium is assumed, such that brine salinity is constrained by temperature (S br = S br (T)).

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 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 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 is the saline density coefficient. Π is the effective permeability, which is a priori depth dependent. g is the acceleration due to gravity, and c, k, and 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, RJW is a dimensionless convection strength tuning parameter and is the density of brine. Outside of the convecting layer, w is zero. 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 We call this scheme "GN13." Similar to RJW14, there are two tuning parameters: 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 where 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."

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(Butler et al., , 2017Eicken, 1992;Saenz & Arrigo, 2012;Vancoppenolle et al., 2006Vancoppenolle et al., , 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 Equation (6) empirically relates the vertical sea ice temperature gradient, T z , to the change in bulk salinity with time, S t , using a desalination intensity parameter, CW . Gravity drainage is moderated by the liquid fraction, , and is switched off if falls below a critical liquid fraction threshold, 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 10.1029/2019JC015791 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 where = −0.0589 • C·kg·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∕ . We keep 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 CW to be a free tuning parameter. We call this convective reformulation of the Cox and Weeks (1988) parameterization "CWconv."

Diffusive Parameterizations
Another stream of authors (Jeffery et al., 2011;Vancoppenolle et al., 2010) proposed representing gravity drainage with a diffusion equation (Figure 1) 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, 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 (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 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 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 where Π is the local permeability, Δ 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."

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 2 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.

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 2 , 18 cm depth (slab 1); ≈500 cm 2 , 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.

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 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 10.1029/2019JC015791 0 (T,S br ) to be equal to one. This term corrects for the increase in brine conductivity, , as the salinity increases relative to the conductivity at the start of freezing, 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 = 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, The uncertainty on S is calculated by Gaussian propagation of the uncertainty in the measured T (±0.1 • C) and estimated Z 0 (±0.5 Ω) (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.

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 To get brine rhodamine concentration we assume the rhodamine is completely dissolved in brine (Appendix B: justifies this assumption) so that 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).

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 equation (1). Of these, RJW14, GN13, and T13 use the Rayleigh number to diagnose convection whereas CWconv uses . 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  Griewank and Notz (2013) with reference parameter values from Griewank and Notz (2015). c  the slow relaxation term, which is incompatible with the transport of solutes other than salt, and take T to be a desalination strength tuning constant. The required boundary condition at the sea ice/ocean interface 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∕ , where C and 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 where F s = −D e · S br ∕ z and F c = −D e · C br ∕ 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).

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 mod ) and measured (S obs ) bulk salinity for N measured layers.
The best tuning parameter or parameter pair (Table 2) was identified as that which produced the minimum d MA .
There are some large differences between our calculated tuning parameters and those presented in the original papers (Table 2). Our lower value for D 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 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.

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 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.

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 3 Evaluation 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.

Statistics of Simulated Bulk Salinity (S) and Bulk Rhodamine Concentration (C) by Comparison With Slab and Wireharp Observational Values
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 10.1029/2019JC015791 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, 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%.

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.

10.1029/2019JC015791
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 We call this scheme "Krembs." We also added a second, simpler adsorption scheme to the model. In this scheme 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 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.