A 1-d Modelling Study of Arctic Sea-ice Salinity

We use a 1-D model to study how salinity evolves in Arctic sea ice. To do so, we first explore how sea-ice surface melt and flooding can be incorporated into the 1-D thermodynamic Semi-Adaptive Multi-phase Sea-Ice Model (SAMSIM) presented by Griewank and Notz (2013). We introduce flooding and a flushing parametrization which treats sea ice as a hydraulic network of horizontal and vertical fluxes. Forcing SAMSIM with 36 years of ERA-interim atmospheric reanalysis data, we obtain a modelled Arctic sea-ice salinity that agrees well with ice-core measurements. The simulations thus allow us to identify the main drivers of the observed mean salinity profile in Arctic sea ice. Our results show a 1.5–4 g kg −1 decrease of bulk salinity via gravity drainage after ice growth has ceased and before flushing sets in, which hinders approximating bulk salinity from ice thickness beyond the first growth season. In our simulations, salin-ity interannual variability of first-year ice is mostly restricted to the top 20 cm. We find that ice thickness, thermal resis-tivity, freshwater column, and stored energy change by less than 5 % on average when the full salinity parametrization is replaced with a prescribed salinity profile.


Introduction
Sea ice is a multiphase material consisting of salty brine, fresh ice, and gas bubbles and is far from static.Brine moves through the ice and across the ice-ocean interface, transporting dissolved tracers such as salt.The thermal properties of sea ice change along with the phase composition; bubbles form, dissolve, and escape into the atmosphere while chemical and biologic processes occur in the brine.Salt is a core component of sea ice as it, along with temperature, dictates the phase composition of sea ice through the liquidus relationship.It also influences the brine density, the chemical properties, the small-scale sea-ice structure, and the vertical stratification of the underlying ocean via salt transport to the mixed layer.Unfortunately, the salinity of sea-ice is an elusive quantity that is difficult to observe.Many open questions related to the salinity evolution can not be answered due to the limited amount and the isolated nature of ice-core measurements, such as to what extent gravity drainage occurs during ice melt, what causes interannual salinity variability, how first-year ice transforms to multiyear ice, and how bulk salinity is linked to ice thickness.To fill these gaps in our understanding, we here study the salinity evolution of Arctic sea ice and quantify the impact of the salinity evolution on various sea-ice properties using an expanded version of the Semi-Adaptive Multi-phase Sea-Ice Model (SAMSIM) introduced in Griewank and Notz (2013).
To do so, SAMSIM needed to be expanded to model seaice surface melt.The surface of melting sea ice is complex and highly heterogeneous.Meltwater flows horizontally through snow and ice into melt ponds and cracks or percolates vertically through the ice.The properties of melting wet snow differ strongly from those of dry fresh snow, and the ice surface also deteriorates during melt and can form a layer of white deteriorated ice which is visually similar to snow (Eicken et al., 2002).All these processes influence albedo.Due to the large influence ice albedo has on sea-ice evolution, the sea-ice modelling community has produced many albedo and melt pond parametrizations (e.g.Flocco and Feltham, 2007;Pedersen et al., 2009), but otherwise surface melt has received very little attention.All 1-D thermodynamic models since Maykut and Untersteiner (1971) have disregarded the physical structure and high gas fraction of the surface during melt and have treated melting sea ice as freshwater ice with modified thermal properties.
Over the last decade, researchers have begun to parametrize the sea-ice salinity evolution (e.g.Vancoppenolle et al., 2006Vancoppenolle et al., , 2007Vancoppenolle et al., , 2009;;Wells et al., 2011;Hunke et al., 2011;Rees Jones and Worster, 2013;Turner et al., 2013) to study the biogeochemical and physical processes in and below sea ice (e.g.Vancoppenolle et al., 2010;Tedesco et al., 2010Tedesco et al., , 2012;;Jeffery et al., 2011;Saenz and Arrigo, 2012;Jardon et al., 2013).Despite these developments, the only published sea-ice model with a fully parameterized salinity evolution is the Louvain-la-Neuve (LIM) 1-D model of Vancoppenolle et al. (2007) based on the 1-D thermodynamic model of Bitz and Lipscomb (1999).Accordingly, many possible approaches for modelling surface melt and parametrizing salinity remain unexplored in 1-D sea-ice models.We introduce new schemes to parametrize surface melt, flooding, and flushing within our 1-D sea-ice model SAMSIM, making it capable of simulating the full growth and melt cycle of sea ice, including the salinity evolution.
We force SAMSIM with Arctic reanalysis data to study the desalination processes and the resulting salinity evolution in the Arctic.This is the first general multiyear model study of sea-ice salinity throughout the Arctic.The only previous model study of sea-ice salinity is the study by Vancoppenolle et al. (2007), which focuses on two ice-core sites of land-fast ice from 1999 to 2001.Model studies are necessary because measurement campaigns can only provide brief glimpses of the full salinity evolution, whereas we can easily explore a far greater diversity of conditions over a longer time frame.The simulated salinity profiles are compared to ice-core measurements to evaluate the model performance.
We have decided to limit the study to the Arctic because flooding and the corresponding snow-ice formation play a large role in the Antarctic.As explained in detail in Sect.2.4.5, we treat the flooding parametrizations currently implemented in SAMSIM as ad hoc solutions only suitable for dealing with isolated and sporadic flooding events.Accordingly, we will refrain from studying Antarctic ice until flooding is better understood.
The final topic we address is how parametrizing the salinity affects various sea-ice properties important to climate models.As sea-ice components of climate models are slowly becoming more sophisticated and modellers have begun to treat sea-ice salinity as a variable instead of a prescribed value or profile (e.g.Vancoppenolle et al., 2009;Turner et al., 2013), it remains unclear how much model performance can be improved by fully parametrizing the temporal salinity evolution and how sophisticated the parametrizations should be to balance the improvements against the increase in computational cost and code complexity.
This paper is organized as follows.In Sect. 2 we detail how surface melt, flooding, and flushing are implemented in SAMSIM.The section ends with a description of the three separate salinity approaches used to parametrize salinity in SAMSIM.In Sect. 3 we conduct an idealized melting experiment to study flushing and to determine how sensitive SAMSIM responds to changes of key parameters.In Sect. 4 we study the salinity evolution of 36 years of simulated sea ice forced with ERA-interim reanalysis data taken from throughout the Arctic.The simulations are split into first-year and multiyear ice, which are analyzed separately and compared to ice-core data.Readers who are primarily interested in the geophysical insights gained by our simulations can understand most of this section without reading Sects. 2 and 3.The final section uses the same atmospheric forcing as Sect. 4 to quantify the impact of the various salinity approaches on quantities relevant to climate models in order to evaluate whether climate models would benefit from a fully parametrized temporal salinity evolution in their sea-ice sub models.

Model description
For the purpose of this paper we expand the SAMSIM model which we first described in Griewank and Notz (2013).SAMSIM is a 1-D column model which employs a semiadaptive grid.In this section we will introduce how SAM-SIM treats surface ablation and processes related to surface melting as well as flooding.
We provide a very brief description of the fundamentals of SAMSIM in Sect.2.1; a detailed description including the underlying equations and numerics can be found in Griewank and Notz (2013).Following the brief description of SAM-SIM we address a small modification of the gravity drainage parametrizations originally presented in Griewank and Notz (2013) in Sect.2.2.Section 2.3 addresses how sea ice melts in reality and in SAMSIM.The final additions to SAMSIM are the parametrizations of flushing and flooding introduced in Sect.2.4.In Sect.2.5 we describe the three salinity set-ups used in SAMSIM.
All parametrizations introduced in this section were designed for SAMSIM.As SAMSIM has some unique characteristics, such as a gas volume fraction (see Griewank and Notz, 2013) none of the proposed parametrizations can be applied in precisely the same way to the commonly used models of Semtner (1976), Bitz and Lipscomb (1999), and Winton (2000).The differences between the models are mostly related to specific definitions of the ice-ocean interface, snow-ice interactions, meltwater formation, and tracer advection.We have made sure to include all assumptions from which the various parametrizations were derived so that corresponding parametrizations for other models can be derived.An evaluation of our new parametrizations is given in Sect.4.3.

SAMSIM
Each layer of SAMSIM is defined by the four fundamental variables mass m, absolute salinity S abs , absolute enthalpy H abs , and thickness z.Absolute values are simply the integral over the mass-weighted bulk salinity S bu and enthalpy H .The solid and liquid mass fractions ψ s and ψ l , as well as the solid, liquid, and gas volume fraction φ s , φ l , and φ g , are derived from the fundamental variables.A salt-free snow layer can exist on the ice, which has a variable density that affects the snow thermal conductivity.However, the only process currently implemented in SAMSIM which affects the snow density is rainfall into snow.This occurs when rain falls while snow is present, during which the snow thickness remains unchanged while the rain displaces some of the previous gas fraction and increases the mass of the snow layer.A full description of the snow and ice thermodynamics is included in Griewank and Notz (2013).
In this paper, we refer to a specific layer by an upper right index counting from top to bottom, with the exception of the snow layer which is marked with "snow".For example, m 6 is the mass of the sixth layer from the surface, m 1 is the mass of the top ice layer, and m snow is the mass of the snow layer.
SAMSIM is the only sea-ice model to employ a semiadaptive grid which grows and shrinks in discrete steps of z 0 at the ice-ocean interface (Griewank and Notz, 2013).However, at the ice-atmosphere boundary it is necessary to have a freely adjustable boundary to deal with incremental surface ablation and snow-to-ice conversion.This is addressed by letting the top ice layer thickness vary freely between 1/2 z 0 and 3/2 z 0 .Once the top ice layer grows thicker than 3/2 z 0 it is split into two layers, the lower layer of the two with a thickness of z 0 .Similarly, when the top ice layer shrinks below 1/2 z 0 it is merged together with the second layer.A sketch of how a grid with three top ice layers evolves during melt is shown in Fig. 1.
This semi-adaptive grid differs in a few crucial aspects from those used in other models such as the one introduced by Bitz and Lipscomb (1999).Firstly, the number of layers of the semi-adaptive grid is not constant and changes with ice thickness (Fig. 1 in Griewank and Notz, 2013), while other models use a fixed amount of layers which grow and shrink with ice thickness.Secondly, the ice-ocean boundary is not defined in the SAMSIM grid, as discussed in Griewank and Notz (2013).Thirdly, the thickness of the upper layers remains constant throughout the run, with the exception of the top layer.Fourthly, as the upper layer boundaries only move in steps of z 0 , there is no numerical diffusion in the upper layers which results from the constant thickness adjustments used in other models.
The short-wave radiation properties of the ice are set with a number of parameters which determine how much radiation is absorbed at the ice surface and how much of the radiation penetrates into the ice and is absorbed in the lower layers.These parameters are the albedo, "alb", the fraction of penetrating short-wave radiation, "pen", and the optical thickness of the ice, κ.Various parametrizations have been proposed which define the optical properties based on the surface temperature, ice thickness, and ablation rates.In SAMSIM the gas volume fraction could also be used to parametrize the optical properties because the number of air bubbles has a large impact on the optical properties of the ice (Light et al., 2008).However, because the focus of this paper is on the salinity evolution, we will use constant values of "alb", "pen", and κ for sea ice to remove a source of variability in the model results (values shown in Table 1).

Modified gravity drainage
We have implemented a slight change to the calculation of the Rayleigh number of the layer i which is used in the gravity drainage parametrizations introduced in Griewank and Notz (2013) as The terms that enter the equation are the standard gravity g, the density difference between the brine in layer i and the lowest layer ρ i , the distance from the layer i to the ocean h i , the thermal diffusivity κ, the dynamic viscosity µ, and the permeability term i .In Griewank and Notz (2013) minimal permeability i = min( i , i+1 , . .., n ) was used as a simplification of the harmonic mean.However, Vancoppenolle et al. (2013) demonstrated that using the minimal permeability instead of the harmonic mean leads to substantially different Rayleigh numbers.Accordingly, we replace the minimal permeability with the harmonic mean in the definition of the Rayleigh number.So instead of i we use the bulk permeability ¯ i for a Darcy flow through a stack of layers, which is given by the harmonic mean overall layers from i to the lowest layer n: where i is the permeability and z i is the thickness of the layer i.
Changing the definition of the Rayleigh number requires the free parameters α and R crit which link the amount of brine to be readjusted, leaving each layer br i ↓ to the Rayleigh number, time step dt, and layer thickness z via br To readjust α and R crit , the same procedure is used as that which initially determined the free parameters in Griewank and Notz (2013).The procedure numerically derives values which lead to the best agreement between modelled salinity and the laboratory measurements of Notz (2005).Two separate sets of measurements and the mean of the two sets are used, resulting in the following free parameter pairings: α = 0.000510, R crit = 7.10; α = 0.000681, R crit = 3.23; α = 0.000584, R crit = 4.89.As in Griewank and Notz (2013) we will use the values optimized to fit the mean of the two measurement sets as the default values: α = 0.000584, R crit = 4.89.In Sect.4.3.3 the effect of the parameter uncertainty of α and R crit on the multiyear salinity profile is addressed.
Updating the Rayleigh number definition has a noticeable effect on the modelled salinity evolution of both the complex and simple gravity drainage parametrizations.However, the qualitative conclusions of Griewank and Notz (2013) and this paper are unaffected by the changed definition of the Rayleigh number.That the qualitative results are unaffected by the change in Rayleigh number definition can be seen by comparing this paper to the results of Griewank (2014), which uses the same simulations but the original Rayleigh number definition of Griewank and Notz (2013).

Surface melt
There are two main difficulties which complicate simulating surface melt in a 1-D thermodynamic sea-ice model.The first is the strong spatial heterogeneity of melting sea ice.Although certain aspects such as melt ponds can be parametrized, there is no way to overcome the fact that a 1-D approximation is less valid for melting sea ice than for growing sea ice.The second major difficulty is that many physical processes which occur at the surface during sea-ice melt are poorly understood.This is especially true for processes which occur at the snow-ice boundary and processes which involve capillary forces in snow or ice.
We have decided against separating the 1-D column into a ponded and non-ponded fraction because this is impossible without violating the core assumption of SAMSIM that each layer is horizontally and vertically homogeneous.A possible compromise would be to couple a 1-D column with a melt pond cover to another 1-D column with no pond, which would come with its own issues of how these columns interact with each other.The classic approach is to implement a melt pond and albedo parametrization which is applied evenly to the column surface without taking any horizontal variability into account.However, we have decided to not introduce such an albedo parametrization for two reasons.Firstly, most albedo parametrizations are not suitable for SAMSIM.For example, some parametrizations change the albedo as an empirical function of surface temperature.If the parametrization assumes that the surface layer is salt free, the parametrization will assume that the surface temperature during melt will always be at 0 • C.However, in SAM-SIM the surface temperature varies during melt depending on the salinity of the top ice layer.Other parametrizations rely on the surface melt speed, which is not a variable in SAM-SIM.Instead, SAMSIM has meltwater formation and surface ablation, which are linked but not identical to the definition of surface ablation used by Bitz and Lipscomb (1999).Secondly, slight albedo changes would overshadow the effects of the sea-ice salinity.If the albedo parametrization were fully physically consistent with SAMSIM then this would be acceptable.However, albedo parametrizations mostly rely on empirical measurements, are intended to improve largescale models, and are ill-suited to determine how the albedo would react to a 5 % increase of gas volume fraction or a 0.1 • C increase of temperature in the top ice layer of SAM-SIM.Including an albedo parametrization would result in a large non-physical source of variability which would greatly complicate interpreting the results.Extending SAMSIM by an albedo parametrization that is compatible with SAMSIM physics remains desirable, however, and will be the subject of future work.For now we simply use a constant value for the ice albedo.
From the measurements taken at the Surface Heat Budget of the Arctic Ocean Project (SHEBA) site, Eicken et al. (2002) identified three stages of melt for Arctic multiyear ice.During stage I melt ponds form, fed by the horizontal transport of melting snow.The snow cover still persists and, while most of the meltwater movement is horizontal, some meltwater drains to the bottom of the ice through cracks and flaws in the ice.Stage II begins when the snow cover has completely melted away.During stage II meltwater moves horizontally until it reaches flaws as well as vertically through the ice.In stage III the flaws have enlarged to the point of ice disintegration.Meltwater moves vertically through the ice and horizontally until it reaches cracks and the edge of the ice flows, and convective overturning occurs in the ice close to the ice-ocean interface.
In SAMSIM, surface melt is implemented by separating melt into two separate stages.The first stage is snow melt, in which snow is converted to slush.This process thins the snow layer by transforming a fraction of the snow into slush, which is then added to the top sea-ice layer as described in Sect.2.3.1.The second stage is surface ablation, in which a fraction of the liquid volume of the top ice layer is designated as meltwater as described in Sect.2.3.2.This meltwater is either transported directly into the ocean or flows through the ice and cracks according to the flushing parametrization introduced in Sect.2.4.2.

Snow melt
The physics of snow is very complex.The snow layer in SAMSIM is intended to simulate only the most basic aspects of snow on sea ice.In contrast to the widely used 1-D thermodynamic sea-ice model of Bitz and Lipscomb (1999), which is implemented in both the Los Alamos (CICE) and the Louvain-la-Neuve sea-ice models, snow does not turn directly into meltwater in SAMSIM.Instead, melted snow from the snow surface percolates downward and accumulates on the sea-ice surface, forming a slush layer of depth B as illustrated in Fig. 2.This snow-to-slush conversion in SAM-SIM is based on two core assumptions.The first assumption is that the snow can only retain a maximum liquid mass fraction (ψ l, max ) which is a function of the snow solid mass fraction.The function we use is which we take from the laboratory study of Coleou and Lesaffre (1998).In Fig. 2 the volume fractions are shown instead of the mass fractions because the volume fractions are proportional to the area depicted.
The second core assumption is that when the liquid water content surpasses the retainable amount, the excess water pools at the bottom of the snow layer, forming a layer of slush.At each time step the depth of the slush layer is deter-mined and then the slush layer is added to the top ice layer.Since the slush layer is merged with the top ice layer as soon as it forms, there is never a slush layer present at the beginning of the following time step.As such, the slush layer is not a physical representation of any physical material but instead a means to transform the model definition of snow into the model definition of sea ice.However, as the model definition of sea ice does not limit the liquid fraction, the sea ice can be in a condition which could be referred to as slush.
Two additional assumptions are required to determine the slush depth which is marked as B in Fig. 2: the gas fraction of the slush φ g, melt and the solid fraction of the slush layer and remaining snow layer.We assume that the solid volume fraction equals the solid fraction of the previous time step and that φ g, melt is a constant.In this paper we set φ g, melt to 20 %, which we base on the measured surface sea-ice densities of Eicken et al. (1995).
Following these assumptions, when the liquid volume fraction of the snow layer exceeds φ l, max , the slush depth B is calculated from the snow solid fraction of the last time step (φ snow s ) and the gas content as (5) As a result, the top ice layer grows thicker by B, and mass and enthalpy are transferred according to the composition of the slush layer.To maintain the solid fraction of the last time step, the snow needs to be reduced in thickness by A as illustrated in Fig. 2. In total, the snow-to-slush conversion shrinks the snow layer by A + B, the total snow and ice column shrinks by A, the top ice layer grows by B, and the snow layer retains its density.
To our current knowledge, the approach of converting snow into slush before it can run off as meltwater is unique.Compared to the standard approach, in which snow melts at the top of the snow layer and immediately runs of as meltwater, our approach leads to a slight delay in the onset of flushing.This delay is because our approach requires the whole snow layer to convert to the model definition of sea-ice via slush formation before runoff occurs.
{ Figure 3. Sketch of meltwater formation caused by surface melting as described in Sect.2.3.2.The white, blue, and grey areas represent the solid, liquid, and gas volume fractions of each model layer (φ s , φ l , and φ g ).The meltwater is located in a film which is z melt thick and located below the surface of the top layer.z melt is determined by the amount of latent heat release necessary to balance the energy difference between the atmospheric heat flux to the surface q atmos and the flux from the surface into the top ice layer q 1 .In reality, sea ice has a varying surface height, which causes the meltwater in the slush to flow into melt ponds.In SAMSIM, by the time the snow layer has melted away, the top model layers that were formed by snow-to-slush conversion are predominantly liquid and salt free but also contain the solid fraction of the meltwater-soaked snow.These top ice layers can be interpreted as a spatial average over melt ponds and snow remnants.As a result the snow melt stage of SAMSIM is shorter than the first melt stage of Eicken et al. (2002) because not all of the latent heat which resided in the snow layer before the onset of melt needs to be released before the snow layer disappears in the model.Although the implemented snow-to-slush conversion neglects many of the finer aspects of snow physics, our approach, by having some interaction between the meltwater which forms at the snow surface with the underlying snow, captures snow melt somewhat more realistically than the standard approach of turning snow directly into meltwater.
Two additional processes also convert snow to slush: flooding as introduced in Sect.2.4.5 and meltwater wicking.Wicking occurs when the top ice layer is so liquid that excess brine seeps into the snow.This process is incorporated into the model as introduced in the following subsection.

Surface ablation
Surface ablation in general refers to an ice-thickness decrease at the surface.Surface ablation is by necessity linked to a flux of melted ice away from the ice surface.In SAMSIM, surface ablation occurs when liquid from the top ice layer is removed via flushing.In this subsection we describe how SAMSIM determines how much liquid is available to be removed from { { { solid liquid gas Figure 4. Formation of meltwater in the top ice layer when φ 1 s < φ s, melt as described in Sect.2.3.2.The thickness of the layer of meltwater ( z melt ) is determined by how much the solid fraction has to be raised to equal φ s, melt .The white, blue, and grey areas represent the solid, liquid, and gas volume fractions of each model layer.
the top layer, what the properties of this liquid is, and how this liquid interacts with the snow layer and the top ice layer.
To describe this clearly we must first clarify how meltwater is defined in SAMSIM.The model definition of meltwater is the liquid in the top layer with the ability to leave the top layer.This ability distinguishes meltwater from the rest of the liquid in the top layer.Otherwise meltwater is identical to the remaining liquid in the top layer (i.e.temperature, salinity, density).Meltwater is assumed to be located on the ice surface in the top sea-ice layer as a thin film.The meltwater film is a part of the top layer, and its thickness is z melt as shown in Figs. 3 and 4.
The amount of meltwater which is present in the top layer is a diagnostic variable which is computed at each time step independently of the amount of meltwater in the previous time step.As the amount of meltwater determines the meltwater film thickness, the thickness is also calculated anew at each time step.
Meltwater can leave the top layer via two processes.The first process is via parametrized flushing, which is detailed in Sect.2.4.2 and 2.4.4.Flushing leads to surface ablation because the thickness of the top layer is reduced by the thickness of the meltwater film when the water flushes away.The second process by which meltwater can leave the top layer is via wicking into the snow layer, as explained at the end of Sect.2.3.1.
SAMSIM relies on three assumptions to diagnose meltwater amount and thickness.The first is that ice melted at the surface of the top sea-ice layer instantly turns into meltwater.The second is that if the solid fraction of the top ice layer sinks below a minimal low value, excess brine turns into meltwater.The third is that over time the gas fraction increases until it reaches the value of φ g, melt .The first two assumptions determine how much meltwater is available in the top layer, while the third assumption influences how thick the melt film is.
SAMSIM determines if melting occurs at the ice surface by analyzing the heat fluxes at the surface.As soon as the sur-face temperature surpasses the freezing temperature given by the bulk salinity of the top ice layer, meltwater can form.The amount of meltwater formed is determined by the amount of latent heat release necessary to balance the energy difference between the atmospheric heat flux to the surface and the flux from the surface into the top ice layer (depicted in Fig. 3).This approach is commonly used in sea-ice thermodynamic models (e.g.Bitz and Lipscomb, 1999) but needs to be adapted to incorporate the varying density and gas fraction of SAMSIM.The discretized diffusive heat flux from the ice surface into the top ice layer is The thermal conductivity of the top ice layer k 1 is a linear combination of the liquid and solid phases, while the gas phase is treated as an insulator.The depth of the meltwater film for a given atmospheric energy flux q atmos is then The second way meltwater can form is when the solid fraction of the top ice layer φ 1 s falls below a minimal low value φ s, melt .When this occurs the solid fraction φ 1 s is rearranged by z melt until φ 1 s reaches φ s, melt , as shown in Fig. 4. From volume conservation it follows that This second way of forming meltwater ensures that meltwater forms before the top ice layer is fully liquid.Not shown in the Fig. 4 is that a similar limit exists on the gas fraction which arises from our third assumption that the gas fraction increases to a specific value over time.If the gas fraction exceeds φ g, melt then the top ice layer is compacted to reduce φ 1 g to φ g, melt , which also slightly increases the density of the top layer.φ g, melt is the same parameter which determines the amount of air captured in the slush during snow melt and is set to 0.2 based on density measurements at the surface of Eicken et al. (1995).To our knowledge there are no measurements from which to estimate φ s, melt .As first guess we assume a value of 0.4, which is slightly above the solid fraction assigned to fresh snow in SAMSIM.If meltwater forms primarily due to low solid fractions, the top ice layer will approach the given values of φ g, melt and φ s, melt over time.
If the meltwater forms due to a low solid fraction while snow is present, the meltwater is assumed to wick up into the snow and creates a slush layer which is then added to the top ice layer again.We refer to this as wicking, and it is similar to snow melt (Fig. 2).The difference between wicking and snow melt is that in wicking the amount of water available to form slush is given by the amount of meltwater present in the top ice layer, while in snow melt the amount is given by how far the liquid fraction of snow exceeds the threshold limit.

Salinity parametrizations
There are three known relevant desalination processes in sea ice: gravity drainage, flushing, and flooding (Notz and Worster, 2009).We addressed how gravity drainage is implemented in SAMSIM in our previous publication (Griewank and Notz, 2013).In this subsection we introduce parametrizations for flushing and flooding, making SAMSIM the second published 1-D model capable of capturing the full salinity evolution.The first model capable of capturing the full salinity cycle is the 1-D LIM sea-ice model of Vancoppenolle et al. (2006).
Parametrizing flushing faces the same challenges as modelling surface melting, namely high horizontal heterogeneity, insufficient data, and a lack of theoretical understanding.No quantitative laboratory studies of flushing have been published to this date and, due to sampling issues and challenging conditions, field studies have been limited to studies of dye dispersion and ice-core salinity (Eicken et al., 2002).The understanding of flooding is even poorer and is limited to the analysis of ice cores which contain flooded snow ice.

Flushing
The first and only published flushing parametrization incorporated in a full thermodynamic sea-ice model by Vancoppenolle et al. (2006) assumes that once the ice reaches a certain permeability, a fraction of the meltwater flows downward through the sea ice and into the ocean below.Although this approach neglects many aspects of flushing, it is able to reproduce field measurements of salinity (Vancoppenolle et al., 2007).In this subsection we will introduce two parametrizations.The complex parametrization attempts to model flushing as a physically consistent hydraulic system, and the simple parametrization is a numerically cheap alternative based on the assumption that the liquid fraction increases towards the surface during surface melt.

Complex flushing
It is known from the field observations of Eicken et al. (2002) that much of the brine movement during flushing occurs horizontally in the upper layers.Once the horizontally flowing meltwater reaches a flaw or crack it drains below the sea ice, which can lead to underwater ice formation (Eicken et al., 1995;Notz et al., 2003).These cracks can also be situated below melt ponds as discussed by Polashenski et al. (2012), who refer to them with the term macroscopic holes.The parametrization of Vancoppenolle et al. (2006) has no explicit treatment of horizontal fluxes.Our goal is to design a flushing parametrization which is as physically consistent as possible in a 1-D model and includes horizontal brine fluxes which are highest when close to the ice surface.Additionally the parametrization should have as few free parameters as possible.The resulting parametrization (sketched in Fig. 5) treats sea ice as a hydraulic network in which each model layer has a vertical and a horizontal hydraulic resistance (R v and R h ).The assumptions on which the parametrization is based are as follows: 1. Cracks always exist in the ice.
2. As we have no data from which to deduce the frequency of these cracks, as a zero-order first guess we assume average horizontal distance between these cracks grows linearly with ice thickness.
3. Once brine reaches such a crack it drains away to the ice-ocean interface without interacting with the underlying ice layers.
4. The vertical resistance represents the resistance to brine flowing from the top to the bottom of a layer.The horizontal resistance represents the resistance that brine needs to overcome to reach a crack.
5. Flushing meltwater flows vertically from layer to layer and horizontally to the cracks.The specific amount for each layer is determined by the hydraulic resistances and the hydraulic head.
6.The hydraulic head is assumed to be equal to the freeboard ζ , resulting in a pressure difference of p = ζρg for the brine density ρ and gravitational constant g.
The resulting parametrization has only a single free parameter β which determines the average distance x to the next crack for a given ice thickness h through x = β • h.
The Darcy flow in a porous medium with a hydraulic resistance of R leads to a mass flux f of for the pressure difference p and liquid density ρ.In SAM-SIM, for each layer i the vertical hydraulic resistance is defined by the permeability , which is a function of the layer's liquid fraction φ i l , the brine viscosity µ, the column area A, and the layer thickness z.SAMSIM uses the permeability function of Freitag (1999), which was derived from measurements of vertical flows.We use it here for both horizontal and vertical permeability.This simplification should not adversely affect our results, since the major simplification lies in the underlying assumption that the permeability is only a function of solid fraction.
To define the horizontal hydraulic resistance we take the average distance to the next crack from our assumptions, resulting in In contrast to the vertical flow area A, which is always 1 m 2 in the column model, the horizontal flow area A i h varies with layer thickness as well as with the geometry of the cracks and resulting flow field.We take A i h to be equal to the vertical layer surface with an area of z i • 1 m.The resulting horizontal and vertical brine fluxes (f h and f v as shown in Fig. 5) are then computed from hydraulic head and resistance.The total resistance over multiple layers is calculated as a sum of parallel and serial resistances, the same method used in resistor ladder circuits.To illustrate how the layers interact, refer to the sketch with six layers shown in Fig. 5 as an example.The lowest layer 6 has by definition no hydraulic resistance.The total resistance of the second lowest layer 5 (R 5 total ) is because R 5 v and R 5 h are connected in parallel.The total resistance over layers 4 and 5 (R 4 total ) is the parallel resistance of R 4 v with the serially connected R 4 h and R 5 total , resulting in Generalizing this for all layers results in which is true for any number of layers.The total amount of flushing brine through the whole ladder circuit shown in Fig. 5 is accordingly The total amount of flushing brine can not exceed the amount of meltwater present in the top ice layer.
The calculated vertical fluxes advect salt and heat from layer to layer using the upstream method, while horizontal fluxes transport both salt and heat directly to the lowest model layer, i.e. the ice-ocean interface.As the thermal profile in melting ice is almost uniform and the brine salinity is linked to temperature, the vertical fluxes lead to a smaller desalination than the horizontal fluxes.
Although the top ice layer can accumulate meltwater faster than it can flush it away, a fully liquid top layer in the model is impossible with the complex flushing parametrization.As the top ice layer becomes more and more liquid, the permeability increases and the horizontal hydraulic resistance of the top ice layer decreases, resulting in a strong horizontal flushing in the top ice layer.This strong flushing removes water from the top layer and prevents the top layer from ever becoming fully liquid.

Complex flushing examples
To illustrate the fluxes which result from the complex flushing parametrization, we apply some numbers to a specific example with six layers as shown in Fig. 5.For this simple thought experiment, layers 1 through 5 are identical with the same permeability and 20 cm thick.Accordingly, the vertical and horizontal resistances of each layer are equal to each other: R The ratio of R h and R v is determined by the free parameter β, the layer thickness z, and the total ice thickness h.In our example we chose z = 0.2 m, which results in h = 1 m because we have five layers of ice.Combined, these result in We can now calculate the ratios of the resulting fluxes for a given value of β, which are shown in Table 2.
For the default value of β = 1, 60 % of the flushing brine would penetrate vertically through all five layers of the ice while 40 % of the flushing brine would flow horizontally until falling through cracks and flaws (row (a) of Table 2).The horizontal fluxes are strongest in the top layer and decrease with depth.A lower value of β would favour horizontal fluxes.Reducing β to 0.2 results in only 18 % of the brine flushing vertically through all five layers, while over 50 % flushes horizontally in the three top layers (row (b) of Table 2).
In the previous example all layers have the same permeability.To illustrate how the complex flushing layer reacts if the lower ice is less permeable, we repeat the same scenario with a higher permeability close to the surface.Specifically, let us assume that the top two layers are 20 times more  5 to illustrate the complex flushing parametrization.All fluxes are given in percent of total flushing ( f h + f 5 v = 100).In (a) and (b), all five layers have the same permeability, while in (c) the top two layers are 20 times more permeable.The free parameter β is changed from the default value of 1.0 to 0.2 in (b).permeable than the lower three.This reduces the percentage of brine that flushes vertically through the whole ice layer from 60 to 15 %, while over 80 % leaves the ice in the top two layers horizontally (compare row (c) to (a) in Table 2).Meanwhile, the horizontal flushing in layers three to five is very small.Less than 5 % of the total brine leaves the ice through cracks and flaws in the lower three layers.If the third layer were impermeable, all flushing would occur horizontally through the top two layers.The scenario of higher permeable upper layers is slightly more realistic than the uniform permeability scenario; however, SAMSIM is run with many more layers and a correspondingly detailed vertical permeability profile.Idealized simulations illustrating how the complex flushing interacts with salinity and thermodynamics are discussed in Sect.3.

Simple flushing
We propose a second, numerically cheaper, parametrization which we will refer to as the simple flushing parametrization.In contrast to the complex parametrization, which calculates brine fluxes that affect salinity via advection, the simple parametrization directly modifies the salinity to fulfill a stability criterion.This stability criterion is based on the simple assumption that the liquid fraction is highest in the top ice layer during melt and decreases into the ice.If this were not the case, the ice below the top layer could become fully liquid.Indeed, fully liquid pools inside the ice have, to our knowledge, never been observed, although rotten ice and slush layers seem to be common during the melt period.This stability criterion is only applied when surface melt occurs and has no affect on the rest of the year when solid fraction www.the-cryosphere.net/9/305/2015/The Cryosphere, 9, 305-329, 2015 is high enough to prohibit liquid from running off as meltwater.
The implementation is as follows.At each time step the meltwater which forms in the top ice layer as explained in Sect.2.3.2 is removed.The salinity of the meltwater is higher than the bulk salinity over the total layer because the solid fraction of the ice is salt free.Accordingly, meltwater removal leads to a reduction of the bulk salinity in the top ice layer.Over time this ensures that the top layer becomes less saline than the second layer.Given that the temperature difference between the top layers is small during surface melt, the second, saltier layer will gradually become more liquid than the fresher top layer.
To ensure that our assumption is fulfilled and the liquid fraction is highest in the top layer, SAMSIM checks each time step if φ 1 l > φ 2 l .When this occurs, the salinity of the second layer is simply reduced by a fixed fraction .This increases the solid fraction while raising the temperature.
The same procedure is then applied to the third layer, to ensure that the second layer is not less liquid than the third layer, and after that to the fourth, fifth etc. until a layer is reached which is less liquid.As long as φ i l > φ i+1 l , the salinity of layer i + 1 will be reduced by the factor .For example if φ 1 l < φ 2 l < φ 3 l > φ 4 l < φ 5 l , the salinity of the second and third layer are reduced while the fourth and fifth remain untouched.

Flooding
Flooding can occur when the weight of snow pushes the ice below the ocean surface, causing ocean water to well up and flood the snow.The resulting frozen mix of snow and ocean water, called snow ice, can be identified by various means in ice cores, from which we know that flooding occurs mainly in the Antarctic and contributes up to 25 % of ice production in certain areas (Jeffries et al., 2001;Maksym and Jeffries, 2001).We base our understanding and treatment of flooding on the work of Maksym andJeffries (2000, 2001) and Jeffries et al. (2001).To readers interested in flooding we recommend the PhD thesis of Maksym (2001).
Although at first glance flooding seems to be the same process as flushing but with a reversed pressure gradient, there are a number of additional uncertainties.Field measurements have shown that a negative freeboard does not automatically lead to flooding although the lower the freeboard, the higher the chance of flooding is.Additionally, very little is known about what happens to the flooded brine once it reaches the ice surface.As flooding occurs at the bottom of the snow mantel, direct observations of flooding are extremely difficult to obtain.Snow metamorphism is in itself a complex process, but the interactions between flooding brine and snow are even more complex and little research has been devoted to this specific issue.Brine movement must occur at the ice surface after or during flooding, because otherwise snow-ice salinities would be higher than the measured values.
As for flushing and gravity drainage, we again developed two separate parametrizations for flooding.However, the two flooding parametrizations are rather similar.We will simply refer to the slightly more sophisticated parametrization as the complex parametrization and the simpler one as the simple flooding parametrization.

Complex flooding
The complex parametrization assumes that during flooding ocean water passes through cracks and channels in the ice to flood the snow layer.The flooding ocean water is assumed not to interact with the brine in the sea ice: Maksym and Jeffries (2001) showed that if flooding resulted in an upward brine displacement through the whole ice, the resulting desalination would quickly turn the ice impermeable.Experiments with SAMSIM reached the same conclusion as Maksym and Jeffries (2001) that upward brine displacement would quickly turn the ice impermeable (experiments not shown).Although the complex flushing parametrization consists partially of vertical flows that displace brine, these only seldom cause the ice to become impermeable for three reasons.Firstly, as a layer becomes less permeable the flushing brine is increasingly diverted horizontally.Secondly, the temperature gradients are much smaller in melting ice so that brine advection leads to less desalination.And thirdly, the ice is usually cooled by the atmosphere during flooding which can compensate the latent heat released during desalination.
The flux of ocean water to the surface is calculated as a Darcy flow driven by the negative freeboard and limited by the permeability of the least permeable model layer.Here we assume that the permeability function of Freitag (1999) provides a useful estimation regardless of the detailed pathways that the ocean water takes through the ice.Although this is a simplification, the major uncertainty stems from the uncertainty in permeability itself and the poor physical understanding of flooding.
Our approach of using the ice permeability to regulate the strength of flooding can lead to a large negative freeboard if the ice layer is impermeable.To avoid this a maximum negative freeboard ζ max is defined.If the freeboard sinks below this threshold, the flux of ocean water necessary to raise the freeboard to the threshold is determined and applied.
The ocean water transported to the ice surface forms a slush layer which is immediately added to the top ice layer at each time step.This is the same approach SAMSIM uses to imitate snow melt and meltwater wicking into the snow layer (described in Sects.2.3.1 and 2.3.2).However, given a snow solid volume fraction of approximately 30-40 %, this approach would result in the flooded slush layer having a very high salinity of roughly 20 g kg −1 , which is inconsistent with measurements.To avoid this high salinity, we assume that the ocean water which floods the snow simultaneously wicks upward and dissolves additional snow into the slush which leads to a freshening of the slush.The ratio of dissolved to flooded snow is assumed to be constant and is defined by an additional free parameter δ.
In this paper we use a value of 5 cm for ζ max , which is based on the freeboard measurements analyzed in Maksym and Jeffries (2000), and we use a value of 0.5 for δ as a preliminary best guess.

Simple flooding
The simple parametrization is the complex parametrization stripped of the permeability-dependent flooding speed and without snow dissolving into the slush layer.The simple parametrization is identical to the complex parametrization if the free parameters are set to ζ max = 0 m and δ = 0.This means that as soon as a negative freeboard develops, flooding sets in right away and no snow is dissolved into the forming slush.

Salinity set-ups
In Sect.2.4 we have presented four parametrizations, two for flushing and two for flooding.Together with the two gravity drainage parametrizations introduced in Griewank and Notz (2013), SAMSIM now has two complete sets of desalination processes.The first set consists of the complex flushing, the complex flooding, and the complex gravity drainage parametrization.The second set of parametrizations consists of the simple flushing, the simple flooding, and the simple gravity drainage parametrization.The parametrizations of the first set all compute brine fluxes which result in salt and heat advection.Accordingly, the rate of salinity change is determined by the strength of brine flow and the salinity gradients between layers.In contrast, the parametrizations of the second set directly adjust the salinity profile to fulfill defined stability criteria.
We will refer to the first set of parametrizations as the complex salinity approach because it consists of the more sophisticated parametrizations which were designed to be as close to reality as possible.The second set will be referred to as the simple approach because the parametrizations included were developed as simpler alternatives to the parametrizations of the complex approach.
The third and final salinity approach employed in this paper prescribes a depth-dependent salinity profile completely independent of the ice properties.The profile used is a crude approximation of measured multiyear ice salinity and is the same profile introduced and used in Griewank and Notz (2013).The profile consists of a linear decrease in bulk salinity from 34 g kg −1 at the ice-ocean interface to 4 g kg −1 at 15 cm above the bottom and a second linear decrease from the 4 g kg −1 at 15 cm above the ice-ocean interface to 0 g kg −1 at the surface.This approach is referred to as the prescribed approach.The prescribed profile is by choice highly idealized so that the prescribed approach provides a stark contrast to the simple and complex approaches.A more realistic profile could have been derived from simulations using the complex approach but we prefer the idealized profile because it is independent of both SAMSIM and the chosen forcing.
An important aspect of the complex parametrization set is that the simulated brine fluxes result in heat fluxes both in the ice and into the ocean.This is most relevant during growth when gravity drainage continually moves colder brine to the ocean while taking up relatively warm ocean water, resulting in a small but steady increase of oceanic heat flux in our limited model domain.Because flushing mostly occurs in ice close to the freezing temperature, the energy lost due to flushing is small.However, these heat fluxes caused by brine flux lead to the complex approach having a different oceanic heat flux than the prescribed and simple approaches.To avoid this change in oceanic heat input when comparing the three salinity approaches against each other, the heat fluxes resulting from gravity drainage and flushing are subtracted from the lowest layer at each time step for the complex approach.This heat flux modification was already applied in Griewank and Notz (2013) to ensure that the various approaches can be compared to each other.

Idealized flushing experiments
In this section we take a closer look at the complex flushing parametrization to study how it interacts with temperature and salinity as well as how sensitively it reacts to various parameters.We prefer to use an idealized set-up, rather than a set-up based on field conditions, for two reasons.The first reason is that in the idealized experiment we can remove all feedbacks and processes not related to flushing.The second reason is that the idealized set-up allows us to chose conditions that highlight how the flushing parametrization interacts with the salinity and thermodynamics of the sea ice.As the full parameter space of all model parameters which interact with flushing in some way is too large to be fully explored in a useful way, we focus on the two parameters which have the strongest effect.The first of these two parameters is β, which determines the linear relationship of average horizontal flow distance to ice thickness in the complex flushing parametrization.The second parameter is the layer thickness z 0 .
The idealized experiment begins with a 1 m thick homogeneous slab of ice with a bulk salinity of 5 g kg −1 and a temperature of roughly −10 • C. The ocean below the ice is at 34 g kg −1 and 0 • C. A constant oceanic heat flux of 15 W m −2 is applied to the bottom while a constant heat influx of 380 W m −2 is applied to the surface.After subtracting the outgoing thermal radiation at 0 • C at the surface, the net heat input into the surface is slightly below 70 W m −2 .The heat fluxes were chosen such that the 1 m slab of ice melts over 1 month, which is the same order of magnitude found in reality.The cold initial temperature was chosen as it highwww.the-cryosphere.net/9/305/2015/The Cryosphere, 9, 305-329, 2015 lights the thermodynamic interactions of the flushing brine.
All brine fluxes that occur in the experiment are caused by flushing as gravity drainage is deactivated and no flooding occurs.
We will first make some general observation of how flushing occurs in the idealized experiment in Sect.3.1 before analysing how the flushing parametrization reacts to β and z in Sect.3.2 and 3.3.

General observations
In the idealized experiment the homogeneous sea-ice slab melts away over 1 month (Fig. 6).The constant surface heat input results in a constant rate of surface ablation.As the initial ice temperature of roughly −10 • C is well below the freezing temperature of the underlying 34 g kg −1 , water-bottom growth occurs over the first 3-4 days.This newly formed ice retains the 34 g kg −1 salinity as no gravity drainage is activated in this simulation.
Flushing commences once the ice surface reaches melting temperature after a few days.The resulting desalination is clearly visible in the salinity profile as well as in the temperature profile (Fig. 6).The downward flushing meltwater quickly desalinates the upper ice, which causes a release of latent heat that warms the desalinated ice to 0 • C.However, after roughly 1 week the flushing stops penetrating downward into the ice and no further desalination occurs in the ice.By comparing the salinity and temperature profiles we can see that the kink in the 3 g kg −1 salinity contour occurs when ice layers with zero salinity are below 0 • C. As freshwater ice below 0 • C is a complete solid, it is impermeable and flushing can not penetrate below this level.This occurs because the temperature in the lower and saline ice cools the freshly desalinated ice layers, while the isothermal desalinated ice transports no heat via thermal diffusions.
Until the impermeable upper layers have melted away after half a month, flushing is restricted to the top layers.Once the impermeable layers have melted away, flushing begins to penetrate into the ice again.As the interior of the ice is by now quite close to the freezing temperature, the newly desalinated layers do not refreeze, and after a few days the ice is fully desalinated.
Two noteworthy secondary effects of flushing occur in the idealized experiment.The first is that while flushing reduces the bulk salinity close to the surface, it also leads to an increase of salinity in the lower ice (visible in the 7 g kg −1 contour of Figs. 6 and 7b-d).This is caused by the positive temperature gradient near the ice-ocean interface, which leads to the vertically flushing brine moving from colder to warmer layers.As the brine is saltier in the colder layers due to the liquidus relationship, salt advection leads to a bulk salinity increase in the lowest ice layers.This effect disappears if gravity drainage is activated (Fig. 7a), which explains why this salinity increase due to flushing has not been observed to our knowledge.To determine if flushing could in principle  1).Plot background is grey.The 3 and 7 g kg −1 contour lines are included.lead to such an increase in salinity if gravity drainage is absent, experiments with a multiphase material, in which both phases have a similar density to inhibit convection, would be required.An additional requirement needed to generate these high salinities close to the ice-ocean interface is that the oceanic heat flux is relatively small so that the flushing parametrization has sufficient time to transport salt into the lower layers before they melt away.
The other noteworthy secondary effect of flushing occurs at the ice-ocean interface.As described in Sect.3, the fresh meltwater which drains through flaws and cracks flows into the lowest model layer.This results in a freshening of the lowest model layer (e.g.layer 6 in Fig. 5).If the lowest layer freezes after it has been freshened by flushing meltwater, it results in a thin layer of low-saline ice close to the ice-ocean interface.This effect is visible in the salinity plots of Figs. 6  and 7, where a thin line of low salinity at the ice-ocean boundary is outlined by the 7 g kg −1 contour line from 0.2 to 0.4 months and once again briefly at 0.5 months.Because this thin layer of ice formed from meltwater, it is less saline than the ice above it.Since it is less saline, the thin layer has a higher solid fraction than the ice above at the same temperature.This leads to a thin sheet of solid freshwater ice below mostly liquid salty ice above.As a consequence, once the ice with low salinity (which is visible as an orange line in the temperature plot of Fig. 6) melts away, the ice above it melts away very quickly due to the low solid fraction.While this ice layer with its low salinity is similar to the false bottoms observed below summer ice, false bottoms occur in nature due to contact of fresh meltwater with sub-zero ocean water, which creates a negative oceanic heat flux.In the idealized experiment the oceanic heat flux is steady and positive, and the formation is dependent on the nonoccurrence of gravity drainage (compare to Fig. 7a).6, experiment description can be found in Sect.3, default settings are listed in Table 1).The 3 and 7 g kg −1 contour lines are included.(a) Gravity drainage is included, which is otherwise disabled in the experiment.

(b)
The ratio of horizontal to vertical hydraulic resistance β is 0.2 instead of 1.0.In (c) the ratio of horizontal to vertical hydraulic resistance β is 5 instead of 1.0.In (d) the vertical spatial resolution z 0 is 2 mm instead of 1 cm, and in (e) the vertical spatial resolution z 0 is 5 cm instead of 1 cm.

Free parameter β
In this subsection we examine the importance of the single free parameter of the flushing parametrization β.We have no definitive physical or model limits on the possible value of β.Based on tracer studies of Eicken et al. (2002), we expect horizontal flows to be on the order of metres.Accordingly, we expect β to be in the single digits, and as a working assumption we set 1 as the default value.A value of 1 assumes the average horizontal travel distance to a crack equals the ice thickness, which implies that the cracks are on average roughly 4 times the ice thickness apart.However, the exact relationship of average travel distance to average crack spacing is a function of the geometric organization of the cracks and the 3-D flow path the meltwater follows.To test the parameter sensitivity around the default value of 1 we repeated the simulation with a value of 0.2 to 5. Because a high β increases the horizontal hydraulic resistance the higher β is, higher values of β cause weaker horizontal fluxes and vice versa, as was shown in the thought ex-periment of Sect.2.4.3.In the idealized experiment the low value of β = 0.2 leads to a delayed onset and depth of flushing in contrast to β = 5 (Fig. 7b and c).The higher value of β increases the salinity at the ice-ocean interface, which results from more brine flushing completely through the ice.The results for β = 0.2, 1, and 5 differ only slightly, indicating that the complex flushing parametrization is much more dependent on the thermodynamics of the ice than the specific value of β.From the idealized experiment we conclude that changing β has the anticipated effect and that the parametrization has a low sensitivity to changes of β close to our default value of 1.This low sensitivity is an advantage for us because although we lack the data to derive the optimal value of β, having a non-optimal estimate of β should only impact our results slightly.

Vertical resolution
Changing the vertical resolution influences the complex flushing parametrization in many ways.The thickness of the top layer has an impact on how SAMSIM calculates meltwater formation, the grid spacing influences heat diffusion and tracer advection, and higher resolution allows more vertical variability of layer properties such as permeability.
The default value for z 0 in the model is 1 cm.As for β we repeated the idealized experiment with a value 5 times lower (i.e. 2 mm) and 5 times larger (5 cm).These values encompasses the practical range of values usable in SAMSIM.
In the idealized experiment, changing the resolution has only a minor effect (see Figs. 6b and 7d, e).The simulations with layer thicknesses of 5 cm and 2 mm are remarkably similar despite the higher resolution run using 25 times more layers.As a result, we do not expect the flushing parametrization to respond strongly to slight changes in vertical resolution.

Summary
The complex flushing parametrization responds weakly to changes of the parametrization parameter β and the model resolution z.Changing β has the expected effect, but no theoretical expectations or data are available to determine the optimal value.Accordingly, the chosen default value of 1.0 is uncertain and may be off by 1 order of magnitude.However, given the low sensitivity to β, even a change of magnitude would not qualitatively change our results.The vertical model resolution has little influence on the parametrized flushing beyond the change in underlying numerics.It is possible that the complex parametrization performs most realistically at a specific layer thickness or that the optimal value of β is resolution dependent, but this can not be determined until more precise data are available.
In this section we study how SAMSIM simulates the salinity evolution in the Arctic using the complex salinity approach and compare the model output with ice-core data.
We have decided to limit the study to the Arctic because flooding and the corresponding snow-ice formation play a large role in the Antarctic.As explained in Sect.2.4.5, we treat the flooding parametrizations currently implemented in SAMSIM as ad hoc solutions only suitable for dealing with isolated and sporadic flooding events.Accordingly, we will refrain from studying Antarctic ice until flooding is better understood.
Although a basic understanding of the salinity evolution has existed for many decades, the main processes driving this desalination still pose many unanswered questions.Using a model has the major advantage of being able to track the evolution consistently over long periods of time, while seaice cores can only provide snapshots.Simulating the salinity evolution with SAMSIM is an exercise in reproducing a vaguely known result of poorly understood origin.We aim to understand the impact and interactions of the various processes better while at the same time discovering the limitations of the developed parametrizations or the existence of neglected relevant processes.

Model set-up
To imitate Arctic conditions we use 3-hourly ERA-interim radiative fluxes and precipitation to provide the surface conditions for SAMSIM.Nine simulations, each forced with ERA-interim reanalysis data taken from one of nine locations spread over the Arctic, are run from July 2005 until December 2009.The coordinates of the chosen locations from south to north are: 70 • E; and 90 • N. A simulation period of 4.5 years was chosen because it covers four yearly cycles of growth and melt, which covers the age of most Arctic sea ice (Lietaer et al., 2011).
SAMSIM also requires oceanic boundary conditions in the form of ocean salinity and oceanic heat flux.Due to the scarcity of oceanic heat flux measurements and for simplicity's sake, all runs share the same prescribed yearly heatflux cycle, which is sinusoidal and based loosely on the heat fluxes (Huwald et al., 2005a) derived from the SHEBA measurements.The oceanic heat flux is highest in autumn (14 W m −2 ) and lowest in spring (0 W m −2 ).Similarly, a standard ocean salinity of 34 g kg −1 is used for all runs.The model settings and parameters used are listed in Table 1.
It is important to state that the boundary conditions we use are not necessarily a realistic approximation of the true conditions at the specific locations and time from which we chose the reanalysis data.Not only are the oceanic heat fluxes a strong approximation, the precision of the reanalysis data is limited by the lack of observations in the Arctic.Additionally, the influences of dynamic processes such as frazil formation, lead opening, melt ponds, and ice drift can not be accounted for in the 1-D SAMSIM model.Given the lack of melt pond formation and lead openings SAMSIM will tend to underestimate the amount of melt compared to reality.

Sample output
To give an example of the model output we have included the salinity evolution of one of the nine simulations for all three salinity approaches (Fig. 8).We chose the simulation forced with reanalysis data from 75 and 145 • W as it has the same forcing during the first growth season as the growth season analyzed in Griewank and Notz (2013).Note that due to the modification to the Rayleigh number (see Sect. 2.2) the salinity evolution of the first growth season shown in Fig. 8 is not identical to the simulated salinity shown in Fig. 9 of Griewank and Notz (2013).
In the sample output the first-year ice survives the first melt season and is followed by 3 years of multiyear ice.The yearly cycle in sea-ice thickness is clearly visible, with strong interannual variations in minimum and maximum ice thickness due to interannual variations in the forcing data, such as snowfall.The complex and simple approaches (Fig. 8a and  b) both create a detailed salinity profile which evolves during growth and melt with large differences from year to year.In contrast, the prescribed approach (Fig. 8c) has neither interannual variability nor a seasonal evolution.As noted in Griewank and Notz (2013), the simple parametrization desalinates slightly stronger during growth, but during the melt season the complex approach loses more salt.In contrast, the prescribed salinity profile results in an increase of bulk salinity over the ice column during melt.

Ice-core data
We begin analyzing the SAMSIM salinity evolution by comparing the output against salinity characteristics derived from ice-core measurements.Despite its drawbacks, taking ice cores is by far the most widespread method of measuring seaice salinity.Gough et al. (2012) provide a thorough overview of statistical and physical sampling issues associated with ice-core salinity measurements.Due to the high horizontal heterogeneity of sea ice we will only use means over multiple ice cores.It is to be expected that the core measurements underestimate the salinity near the ocean interface due to brine loss (Notz and Worster, 2008).
After over a century of sporadic measurement campaigns beginning with Nansen's Fram expedition, the observational record of Arctic sea-ice salinity is sparse in time and space and no comprehensive compilation of the conducted measurements has been published in the last decades (e.g.Weeks and Lee, 1958;Cox and Weeks, 1974;Nakawo and Sinha, 1981;Eicken et al., 1995).We do not attempt to provide a rigorous comparison of model versus field data in this paper.Instead, we select three characteristic traits of sea-ice salinity to compare SAMSIM's results against.The three traits we compare against are the link between bulk salinity and ice thickness, the first-year salinity evolution from January to June, and the mean multiyear salinity profile from May to September.

Bulk salinity against thickness
The first trait we selected is the link between salinity and thickness which was studied by Cox and Weeks (1974) and Kovacs (1997).For the single growth season studied in Griewank and Notz (2013) the model results agreed well with the fit of Kovacs (1997) for first-year ice up to 2 m.
We separate first-year from multiyear ice before comparing the bulk salinity against thickness (Fig. 9).One simulation was singled out and highlighted, allowing the reader to track the progress over 4 years as the first-year ice turns into multiyear ice and becomes less saline over time.The simulation which was singled out is the same simulation as shown in Fig. 8.
Both first-year and multiyear ice show a distinctly different behaviour during growth and melt.The gradual transition from growth to melt is visible as a drop in bulk salinity at a constant thickness.A closer examination reveals that a slight thickness increase is visible in many simulations before ablation sets in.This bump in ice thickness arises from SAM-SIM's definition of sea ice, which includes melting snow that has turned into slush (for details see Sect.2.3.1).That this little bump appears at the end of the downward drop signals that until then no flushing has occurred.From that we can conclude that gravity drainage causes the drop in salinity.
Ice thinner than 20 cm has a wide spread in bulk salinity caused by melting and flooding at the onset of the growth season.The simulated first-year ice thicker than 20 cm agrees well with the empirical results of Cox and Weeks (1974) and Kovacs (1997) during growth, with the model having only a slightly higher salinity.This bias is especially high for ice thinner than 0.5 m, which may be partially due to the fact that the underestimation of bulk salinity due to brine loss is higher for thin cores.After the onset of melt the bulk salinities are comparable to the estimates of Cox and Weeks (1974), which were based on a limited amount of cores that were at least 1 m thick.
As expected, multiyear sea ice shows a much smaller range of bulk salinities (Fig. 9b).During growth the bulk salinities show no coherent dependence on thickness, but during melt there appears to be a slight linear dependence on thickness.This is not far off from the estimation of Cox and Weeks (1974).
Both the modelled first-year and multiyear profiles are almost completely salt free at the end of the melt season.Neither Cox and Weeks (1974) or Kovacs (1997) included ice cores of such thin ice during melt, so we can not conclude from our comparison if this model behaviour agrees with reality.However, it is plausible that the 1-D nature of SAM-SIM, which is built on the assumption that ice layers are totally homogeneous and all brine pockets are connected, would lead to an overestimation of desalination during flushing.
In conclusion, the modelled thickness-salinity relationship of growing first-year ice agrees well with the empirical fits to measurements of both Cox and Weeks (1974) and Kovacs (1997).For growing multiyear ice there is no oneto-one relationship between thickness and salinity, though growing multiyear ice tends to to be less salty the thicker it gets.The transition from growing to melting ice leads to a loss in bulk salinity at a constant thickness which is caused by  Kovacs (1997) for ice up to 2 m.The red dashed lines mark the empirical linear relationships found by Cox and Weeks (1974) for growing (upper lines) and melting Arctic ice (lower line).
gravity drainage in the warming ice.Both melting first-year and multiyear ice show a weak linear dependence of salinity on thickness.In our simulations, the ice loses almost all its salt during melting; hence its mean salinity after the re-onset of growth is strongly affected by the salinity evolution of the newly forming ice.

First-year salinity evolution
The second trait of the modelled salinity we evaluate with core data is the evolution of first-year ice salinity from January until June.A longer time frame was not possible due to data availability; the period nonetheless allows us to study the salinity changes after gravity drainage is mostly restricted to the lower layers.We use the ice-core data taken as part of the Seasonal Ice Zone Observing Network and the Alaska Ocean Observing System by the sea-ice research group at the Geophysical Institute at the University of Fairbanks from 1999 to 2011 (Eicken et al., 2012).The great advantage of these measurements, other than the sheer number of cores taken, is that by measuring repeatedly over a decade a large spread of conditions were captured.After rejecting all cores which did not include an ice thickness measurement or contained gaps in the salinity profile, a total of 86 first-year profiles remained between January and June.
The comparison of the model salinity with the Barrow cores is not ideal because SAMSIM is forced with conditions from throughout the Arctic, while the cores were all taken close to the Alaskan coast as part of an ongoing effort to understand and alleviate the impact of changing sea-ice on the human settlements along the coast (Druckenmiller et al., 2009).Ideally we would force our model with the forcing experienced by the ice measured at Barrow.This is not possible for a number of reasons.Firstly, we have no measurements of the oceanic heat flux.Secondly, although the ice was measured in Barrow we do not know where it was before the core was extracted.Most of the ice will likely have formed near the extraction points, but as illustrated by the multiyear ice cores taken in a region which is ice free in summer, there is a substantial amount of drift.Thirdly, we do not know when the ice was formed.The ice could have been formed during the initial freeze-up in fall, or later on in a lead or polynya.And lastly, due to the uncertainty in reanalysis data and the high variability in snow depth, we could not be certain that applying reanalysis forcing taken from the exact point where the cores grew would be correct.However, we do have a number of reasons to believe that the model-data comparison is useful.Firstly, the cores are taken over 12 years.This means that interannual variability will ensure that ice grown under a range of conditions was measured.Secondly, we show in the subsection on interannual salinity variability that the salinity variations resulting from atmospheric conditions are strongest in the uppermost 20 cm (Sect.4.5).Because of this, we believe that the comparison should work well for the lower 80 % of the ice.
To compare the core profiles against the model profiles, both are first normalized to a depth of 0 to 1 before averaging over time.Often the salinity measurements did not extend all the way to the bottom of the ice, in which case the lowest measurement was extrapolated downwards.This extrapolation will contribute to the underestimation of salinity at the ice-ocean interface common to ice cores.We group the 86 core measurements into three bins of similar size based on the dates they were taken.The first bin spans from January to March (27 cores), the second from April to May (29 cores), and the final bin contains the remaining 29 cores taken in June.
As expected, even though the core profiles have a sharp increase of salinity at the ice-ocean interface they are still less saline at the ice-ocean boundary than SAMSIM (Fig. 10).Other than the top and bottom 10 % of the ice thickness, the simulated salinity profiles and the Barrow cores never differ by more than 2 g kg −1 , which is in itself a mentionable model feat.and first-year ice from reanalysis forced simulations using the complex brine dynamic parametrizations (b).Both were averaged from January to March (1-3), April to May (4-5), and over June (6).
Other than the general agreement, this comparison highlights some limitations of SAMSIM's complex salinity approach.One of these limitations is that flushing and snow melt by design lead to a zero salinity at the surface once surface melt commences.Accordingly, the June SAMSIM profile is completely salt free at the surface while the core data show a salinity of roughly 1 g kg −1 at the surface (Fig. 10).
This total desalination at the surface is rooted in two of SAMSIM's design choices.The first design choice is that the snow layer in SAMSIM has zero salinity and that melting snow forms slush which is treated as sea ice.Accordingly, when snow melts the top ice layer will consist of melted snow slush and be absolutely salt free.The second design choice which leads to zero salinity at the surface is the implementation of flushing in SAMSIM.One of the core assumptions of the complex flushing parametrization is that the meltwater leaving the top ice layer has a brine salinity determined by the liquidus relationship.Accordingly, as the brine salinity of the top ice layer is by definition always higher than the bulk salinity of the top ice layer, flushing always results in a salinity decrease at the surface.This desalination quickly desalinates the surface once flushing commences as shown by the idealized flushing experiments (Sect.3).While the freshly desalinated ice can freeze solid and thus inhibit any further flushing, this can only occur in the ice if the underlying ice is sufficiently cold as in the idealized example.At the surface this could also occur but only if there were a negative atmospheric heat flux to remove sufficient energy from the top layer to overcome the latent heat released during freezing.
The second distinct difference between model and core salinity is that SAMSIM has a high surface salinity with a very strong salinity gradient before the onset of melt (profiles from January to May, Fig. 10b).The sharp salinity gradient which occurs in the top few model layers could be a numerical artifact arising from SAMSIM's semi-adaptive grid.No matter which resolution is used, the initial ice growth occurs when only a few layers are active.This issue was investigated in our previous paper when comparing to freezing plate experiments conducted in the lab, but available data were insufficient to make any conclusions (Griewank and Notz, 2013).A different explanation is snow wicking, a process which transfers some of the surface salinity into the snow layer.In the model wicking only occurs when meltwater forms in the top ice layer beneath snow.The discrepancy between model and data at the surface could also arise from the neglect of frazil or pancake ice formation in SAMSIM.In frazil and pancake ice the wave motion and turbulence cause brine motion not captured by the gravity drainage parametrization which could desalinate the initial ice before it freezes into a static structure.
The third discrepancy between the cores and SAMSIM is that the bulk salinity in the upper 40 % does not change substantially from the period of January-March to that of April-May in the model.There are many possible explanations for this discrepancy, such as the non-ideal comparison itself (see second paragraph Sect.4.3.2),insufficient simulations or core measurements, and errors of the core-salinity measurements.Another explanation is that the model is unable to simulate the salinity evolution correctly close to the surface during winter.A likely candidate to explain that the salinity remains constant near the surface is that the gravity drainage parametrization desalinates too quickly during growth.The modelled salinity is quickly reduced to 5 g kg −1 after which it stabilizes, instead of a weaker initial desalination followed by a gradual desalination over time (Fig. 10).The neglect of frazil and pancake ice formation in SAMSIM could again be an issue since turbulent conditions during the initial freeze-up would influence both the microstructure and permeability of the surface ice.It is also possible that the freeboard plays an important role, and that brine from above the waterline drains away by an unknown mixture of gravity drainage or flushing.The differences between the cores and SAMSIM as well as our poor understanding of what happens during flooding indicate that unknown, yet relevant, brine movements may occur at the ice-snow interface.It is also possible that the gravity drainage parametrization has some limitations.Despite the indirect model-to-data comparison and the three discrepancies in the salinity evolution discussed, SAMSIM successfully captures the general shape and magnitude of the three core-derived salinity profiles.

Multiyear salinity profile
The final and best-documented trait we select to compare is the mean multiyear salinity profile.The most widely used multiyear profile in the sea-ice modelling community is www.the-cryosphere.net/9/305/2015/The Cryosphere, 9, 305-329, 2015 May to September mean of vertically normalized multiyear salinity profiles of reanalysis forced simulations using the complex brine dynamic parametrizations.Schwarzacher 59 refers to the fitted profile of Schwarzacher (1959), and Barrow cores refers to the multiyear ice cores taken by the Alaska Ocean Observing System from 1999 to 2011 (Eicken et al., 2012).The R crit , α spread shows the SAMSIM profile using the two non-default values of the gravity drainage parameters obtained from the optimization process (see Sect. 2.2).Left line: α = 0.000681, R crit = 3.23.Right line: α = 0.000510, R crit = 7.10.The area between the two simulations is shaded in light grey.
based on 40 ice cores taken at the drifting ice station A in 1958 (Schwarzacher, 1959) from May to September.Although later studies have incorporated additional measurements (e.g.Cox and Weeks, 1974;Eicken et al., 1995), the basic shape has remained similar.The fitted bulk salinity profile of Schwarzacher (1959) on a normalized vertical coordinate z from zero to one, is used in the 1-D models of Maykut and Untersteiner (1971) and Bitz and Lipscomb (1999).As the Schwarzacher cores were all taken from May to September, and the eight multiyear cores from Barrow were also taken in summer, we compare the normalized Barrow cores and Schwarzacher profile against the mean of SAMSIM from May to September.We did not compare directly to the Schwarzacher data as they were not easily available and the data displayed in Schwarzacher (1959) were not regularized before averaging.
Although the fitted Schwarzacher profile has a 3.2 g kg −1 salinity at the ice-ocean interface (Fig. 11), an increase is clearly visible in the measurements, similar to the salinity increase of the eight multiyear salinity cores taken at Barrow.Due to this ignored increase and the repeatedly mentioned salinity loss in cores, we only compare to the upper 90 % the Schwarzacher profile.Although this comparison of SAMSIM to field data is far from perfect, it is the closest we can come to evaluating the flushing parametrization until controlled laboratory measurements are available.
We compare the May-to-September mean of all normalized multiyear SAMSIM profiles to the profile of Schwarzacher (1959) and the Barrow cores, which all share a similar magnitude and shape in the upper 90 % (Fig. 11).Both SAMSIM and the Barrow cores have a slight maximum at a depth of 40 %, which indicates that the complex flushing parametrization predicts the desalination depth reasonably correctly.The good agreement between SAMSIM and ice-core data is a very positive result given that the complex flushing parametrization contains large parameter uncertainties and was developed from scratch without any data available to tune the free parameter β.
Between the depth fractions of 0.5 and 0.8, SAMSIM and the Barrow cores show a slight salinity decrease with depth while the Schwarzacher profile has a slight increase (Fig. 11).The differences between the model and the ice-core data are of similar magnitude to the differences caused by different values of α and R crit obtained from the optimization process mentioned in Sect.2.2.The Barrow cores and SAMSIM both have a sharp salinity increase in the lowest 10 %.That the model is saltier at the ice-ocean boundary is expected due to brine loss during coring and a lower spatial resolution of the measurements compared to the model.

Summary
According to SAMSIM there is a clear link between ice thickness and bulk salinity in growing first-year ice as described by Kovacs (1997).However, after the ice stops growing, gravity drainage in the warming ice causes a thickness independent desalination.Both melting first-year and multiyear ice show an approximately linear dependence of bulk salinity on ice thickness as suggested by Cox and Weeks (1974).The modelled ice loses almost all its salinity, a feature against which we do not have any core data to evaluate.The mean multiyear salinity profile of SAMSIM from May to September agrees well with the core data of Schwarzacher (1959) and from Barrow.The salinity evolution in first-year ice in SAMSIM is comparable to ice-core measurements at Barrow (Eicken et al., 2012).However, in contrast to the Barrow core data, the modelled salinity close to the ice surface remains constant from the period of January-March to that of April-May, indicating that in reality, brine fluxes occur close to the surface and are poorly captured by the complex set of parametrizations.
All comparisons between SAMSIM and ice-core data show that SAMSIM captures the general salinity evolution well, both qualitatively and quantitatively.Keep in mind that no tuning was used to reach these results and that all parametrizations were developed without any field data.Additionally, all parametrizations were developed separately, with no regard to possible interactions.
So far we have only evaluated characteristics of the Arctic simulations that we could compare against ice cores.From the comparison to ice cores we conclude that our parametrizations and understanding of desalination processes are sufficient to use SAMSIM as a tool to study Arctic sea ice beyond reproducing ice-core salinity.

Mean salinity evolution
In this subsection we analyze the mean salinity evolution of the complex approach.In total, the model simulations yield 36 years of sea-ice growth and melt.Of those 36 years, 21 years are multiyear ice and 15 are first-year ice.Of the 15 years of first-year ice, 8 years end in open water while 7 form multiyear ice in the following year.
To process and visualize the salinity evolution we first normalize the depth of all salinity profiles of the model output between 0 and 1.This allows averaging over multiple normalized profiles and it simplifies comparing profiles of varying thicknesses.To resolve the mean annual cycle we sort all first-year and multiyear profiles into monthly bins beginning in September, which we then average (Fig. 12).A side effect of this averaging approach is that when there is no ice in the model output, this output does not affect the mean salinity profile.As a consequence, the mean August profile consists mostly of first-year ice which will turn into multiyear ice the following year, and there is a smooth transition from the August first-year profile to the September multiyear profile.This selection effect is clearly visible when comparing the mean ice thickness of all first-year simulations excluding ice-free output against the mean thickness of first-year ice which turns into multiyear ice next September (Fig. 13).
During the growth season the salinity of the first-year ice decreases to 5 g kg −1 after about 2 months with a sharp increase to 10 g kg −1 in the upper 5 % of the ice thickness (Fig. 12a).The salinity profile remains pretty stable between November and April, followed by a slight desalination in May at the onset of melt.The desalination accelerates during June and July until the upper 80 % of the ice has a very low salinity below 2 g kg −1 (Fig. 12b).The influence of flushing is clearly visible in the almost total loss of salt at the surface from June onwards.Although there is only little and indirect experimental evidence of gravity drainage occurring as the ice warms (e.g.Widell et al., 2006;Jardon et al., 2013) the salinity reduction in the lower half of the ice from April to June shows that gravity drainage is active in SAMSIM during the onset of melt.This desalination is consistent with results from idealized experiments we conducted that show a reduction of bulk salinity from above 5 to below 3 g kg −1 from gravity drainage when sea-ice begins to warm (Griewank and Notz, 2013).
At the end of the melt season the multiyear ice salinity is lowest (Fig. 12c).While the surface salinity remains low the newly formed ice at the bottom retains over 5 g kg −1 .During the melt season the lower half of the ice is desalinated by gravity drainage while flushing maintains the low surface salinity.That this desalination is not only due to the loss of the saltier lower layers through melt is visible in the curve that develops in the lower half of the normalized profile as flushing by itself would lead to an increase in salin- Figure 14.Yearly mean first-year (fy) and multiyear (my) sea-ice salinity profiles of SAMSIM using the complex parametrization.The fitted analytical functions of the profiles listed in Sect.4.4 are added in orange.Although the profile of Schwarzacher (1959) is summer biased (see Sect. 4.3.3),we have included it as a reference.ity (Fig. 6).That gravity drainage can act in such a manner is visible in the idealized experiment in which gravity drainage was enabled (Fig. 7a).This curve is also visible in the Barrow core data shown in Fig. 11.With the exception of the gravity drainage during melt, the overall multiyear salinity agrees well with expectations already voiced by Cox and Weeks (1974).
For readers interested in analytical approximations of the mean first-year and multiyear profiles we offer two functions, S bu, fy (z) and S bu, my (z).Both are a function of the normalized ice depth 0 ≤ z ≤ 1 and are shown in Fig. 14  The transition from first-year to multiyear ice over the melt season can be approximated by a time-dependent combination of the two profiles in the form where t = 0 at the beginning of the melt season in June and t = 1 at the onset of growth in September.

Variability
While the previous subsection studied the mean salinity properties, in this subsection we will take a brief look at the salinity variability in SAMSIM using the complex approach.The model variability arises from two sources, the main one being the atmospheric forcing.Although the location at which the reanalysis data was selected has the largest impact, interannual variability ensures that all 36 years of simulated sea ice have a unique forcing.The second source for variability is the initial ice conditions at the beginning of the growth season.This second source only applies to the 21 years of multiyear ice since all first-year ice grows from ice-free water.The variance of the model can not be directly compared to ice-core variability, because the variability in ice cores additionally contains a large amount of variability due to small-scale horizontal heterogeneity (Gough et al., 2012).
To visualize the variability we have plotted all normalized salinity profiles at two dates in time as well as the mean overall profiles at that time point in Fig. 15.With few exceptions the first-year ice only deviates a few g kg −1 from the mean in the lowest 80 % of the ice.However, at the surface the spread is much higher, with values reaching from 0 to above 10 g kg −1 (see Fig. 15a and b).There are two main reasons for the higher variability at the surface.The first is that after 10-20 cm of ice has formed, the variability of the atmospheric forcing is severely dampened before it reaches the ice-ocean interface.As a result, the ice formed after the initial 10-20 cm grows under roughly similar conditions in all simulations.The second reason is that flooding and flushing both occur mainly at the surface of the ice.That such a similar high variability near the surface is not visible in the multiyear ice is because both processes are far less likely to occur in multiyear ice during the winter than in first-year ice.Farther south, where first-year ice seldom survives the melt season, rainfall and above-freezing surface temperatures occur during the growth season, both of which can cause flushing.As the first-year ice is less thick, strong snowfall that slows ice growth can lead to flooding more easily than in multiyear ice.
As all multiyear ice has experienced at least one melt season, it is not surprising that multiyear simulations have a salinity of zero at the surface (Fig. 15c and d).That all 21 years have zero surface salinity shows that flooding of multiyear ice does not occur in any of the simulations.Most of the variability in multiyear ice arises from the different ice thickness and salinity of the ice at the end of the melt season.The sudden salinity increases with depth arise from sudden quick growth in the beginning of the growth season beneath almost completely desalinated ice for many simulations between 0.2 and 0.6 in the November profiles (Fig. 15c).This growth can be quicker than in first-year ice of similar thickness due to the following reasons.The first reason is that by the time first-year ice reaches the same thickness, it has likely accumulated an insulating snow layer which slows ice growth.Secondly, the fresher multiyear ice has a higher thermal conductivity and lower thermal capacity which enhances heat transport from the ice-ocean interface to the iceatmosphere boundary.
Over the next half-year the profiles are smoothed out and the salinity sinks to 7 g kg −1 or lower except in the lowest 10 % (Fig. 15d).Visible in both first-year ice and multiyear ice is that the salinity in the lowest layers is higher in November during ice growth than in April.
In conclusion, the variability in first-year ice is strongest at the surface and arises from the atmospheric forcing, while the variability in multiyear ice is mostly due to the thickness of the ice at the beginning of the growth season.A third possible source of variance is the variation in the oceanic heat flux.This is not included in this study as all simulations share the same prescribed annual cycle of oceanic heat flux.

Impact of parametrizing salinity
While the previous section focused on the salinity evolution and the processes which drive it, this section aims to quantify how parametrizing salinity affects sea-ice properties relevant to the climate system.We address this question, which is highly relevant to modellers seeking to improve climate models, by using the same runs used in the previous section (see Sect. 4.1).In this paper we only study the physical properties of sea-ice.Biogeochemical properties and feedbacks can not be assessed with SAMSIM currently.
To asses the total impact of parametrizing salinity in a climate model it is not sufficient to quantify the impact on the sea ice itself.It is also necessary to determine resulting feedbacks with the ocean and atmosphere.So far the only coupled model featuring a partially parametrized salinity is the NEMO-LIM model, which uses a prescribed atmospheric forcing.Using the NEMO-LIM model, Vancoppenolle et al. (2009) found that the large-scale sea-ice mass balance and the upper-ocean characteristics are quite sensitive to sea-ice salinity.Salinity variations introduced to NEMO-LIM increased sea ice volume by up to 28 % in the Southern Hemisphere because changes to ice-ocean interactions stabilized the ocean, leading to a reduced oceanic heat flux.In the Arctic the ocean stratification was not influenced by the implemented sea-ice variations; however, Vancoppenolle et al. (2009) discovered increases in ice thickness of up to 1 m due to changes of the sea-ice thermal properties.
From Vancoppenolle et al. (2009) we conclude that in the Arctic the oceanic feedbacks will be small due to the stable stratification of the Arctic Ocean.Although the atmospheric feedbacks remain unknown, we can use SAMSIM's more advanced salinity parametrizations with a much higher spatial and temporal resolution to take a more detailed look at how the salinity evolution affects the sea ice.
To quantify the impact of parametrizing salinity we compare quantities of the nine reanalysis forced simulations using the three salinity approaches introduced in Sect.2.5.The specific quantities we use based on their importance for the climate system are the same four used in Griewank and Notz (2013).These are the ice thickness, the freshwater column stored in the ice and snow, the thermal resistance R th , and the total enthalpy H integrated over the whole ice and snow column.Each of the nine runs is evaluated separately over the full 4.5 simulation years to ensure that opposing biases at different locations do not average out.
The metrics we use to compare the time-dependent quantities against each other are a time-integrated ratio and a timeintegrated, weighted absolute difference.The ratio r of the quantity x i (t) using the salinity approach i to the same quantity using the different salinity approach x j (t) over the simulated 4.5 years is calculated as r = t=4.5 a t=0 x i (t) dt t=4.5 a t=0 x j (t) dt .
(20) Figure 16.Ratios of the time-integrated ice thickness, freshwater column, thermal resistance R th , and enthalpy H of the simple and the prescribed SAMSIM salinity approach compared to the complex approach (details in Sect.5).The ratios were calculated separately for each of the nine reanalysis forced simulations over 4.5 years.
Each dot shows the ratio of a specific simulation, while the lines show the mean overall runs and quantities.
The second metric used, the weighted absolute difference "d", is determined by d = t=4.5 a t=0 x i (t) − x j (t) dt t=4.5 a t=0 x j (t) dt ( 21) and is a measure of how large the differences are between the two quantities at each time step compared to the total value of the second quantity.The ratio is chosen to indicate if and by how much x i is greater or smaller than x j over time, while the absolute difference is chosen to detect compensating errors not apparent in the ratio.We quantify the impact by comparing the simple and prescribed approach against the complex approach.
The computed ratios for each simulation reveal that the prescribed approach with few exceptions leads to a lower ice thickness, freshwater column, thermal resistance, and total enthalpy than the complex approach (Fig. 16).Ratios range from 0.90 to 1.05.The mean of all ratios and quantities of the prescribed approach is 0.975; accordingly, the quantities of the complex approach are 2.5 % higher on average.The ratios of the simple approach have a slightly lower spread and are on average higher with a mean of 1.012.So on average the simple approach overestimates half as much as the prescribed approach underestimates.
The absolute differences paint a similar picture, with the prescribed approach having a slightly larger spread with differences up to 12 % (Fig. 17).On average the simple approach has lower differences with a mean of 3.3 % in comparison to the prescribed mean of 4.5 %.Because the absolute differences are roughly 2 times larger than the ratios, we can deduce that roughly half of the discrepancy between two simulations stems from a bias in one direction.
Given that the prescribed approach does not distinguish growing from melting ice and that the prescribed profile was not optimized or tuned in any way, the simulated ice properties using the prescribed approach are unexpectedly close to the complex approach.We also expected the prescribed ap- .Time-integrated absolute differences of the ice thickness, freshwater column, thermal resistance R th , and enthalpy H of simulations using the simple and prescribed SAMSIM salinity approach compared to simulations using the complex approach (details in Sect.5).The absolute differences were calculated separately for each of the nine reanalysis forced simulations over 4.5 years.
Each dot shows the ratio of a specific simulation, while the lines show the mean overall runs and quantities.
proach to have a wider spread when compared to the complex approach, because the prescribed approach treats all ice the same regardless of its history while the complex approach is dependent on previous conditions (as visible in Fig. 8).

Summary and conclusions
We have incorporated surface melt, flooding, and flushing into SAMSIM.In contrast to the thermodynamic models derived from Maykut and Untersteiner (1971), such as Bitz and Lipscomb (1999) and Huwald et al. (2005b), surface melt in SAMSIM is implemented as a two-stage process.The first stage is the conversion of snow to slush followed by the second stage of surface ablation by meltwater runoff.All desalination processes are parametrized in two different ways in SAMSIM.The complex parametrizations calculate brine fluxes and are physically consistent, while the simple parametrizations attempt to imitate the effects of the complex parametrizations with less numerical overhead.SAMSIM is the only 1-D thermodynamic sea-ice model other than the 1-D LIM model of Vancoppenolle et al. (2007) which has a fully prognostic salinity.In contrast to the flushing parametrization of Vancoppenolle et al. (2007), the complex flushing parametrization of SAMSIM explicitly includes both horizontal and vertical brine movements.A detailed discussion of why the complex gravity drainage parametrization of SAMSIM agrees better than the gravity drainage of LIM 1-D with both theoretical and numerical expectations is included in Griewank and Notz (2013).The complex flooding parametrization based on the results of Maksym and Jeffries (2000) is an ad hoc solution as the current understanding of flooding is insufficient to develop a more realistic parametrization.Nevertheless, SAMSIM is the first 1-D model to include flooding as well as flushing and gravity drainage, and the flooding parametrization does capture the basics of flooding and produces snow ice with reasonable salinities in a physically consistent manner.
Under idealized conditions, the complex flushing parametrization leads to an increase of salinity close to the ice-ocean interface if gravity drainage is deactivated.Although we do not have data available to determine optimal values of the ratio of vertical to horizontal hydraulic resistance β, our idealized experiments show that the flushing parametrization is only weakly sensitive to changes close to the default values.The vertical resolution of SAMSIM also only has a small impact on the flushing parametrization.
We study the salinity evolution of Arctic sea ice using 36 years of SAMSIM output.To imitate Arctic conditions we force SAMSIM with ERA-interim reanalysis precipitation and radiation fluxes from throughout the Arctic.The 36 years are separated into 15 years of first-year and 21 years of multiyear sea-ice and then compared against ice-core data.The mean multiyear salinity profile of Schwarzacher (1959) and the salinity evolution of first-year ice cores from Barrow, Alaska, agree well with SAMSIM simulations.However, while the first-year ice-core salinity at the surface decreases from January to May, the modelled salinity at the surface remains constant until the onset of melt.This discrepancy indicates that brine fluxes close to the ice-snow boundary are captured poorly by SAMSIM.Possible reasons for this discrepancy are discussed in detail in Sect.4.3.2.
We deduce from the 36 years of simulated sea-ice that ice thickness is a good indicator of bulk salinity for growing first-year ice.The model results agree well with the empirical results of Cox and Weeks (1974) and Kovacs (1997).That the modelled bulk salinities of thin ice are higher than the icecore data is at least partially due to the fact that brine loss during coring is especially high from thin and more saline ice.The transition from growth to melt is accompanied by a 1.5-4 g kg −1 reduction of bulk salinity caused by gravity drainage before the onset of flushing.This onset of gravity drainage as the ice warms is consistent with earlier findings by Griewank and Notz (2013) and Jardon et al. (2013).The onset contradicts the general melt evolution depicted by Eicken et al. (2002) in which gravity drainage sets in at the end of the melt season.In general, thicker multiyear ice tends to be fresher, but during growth the bulk salinity increases with thickness.During melt both multiyear and first-year ice have a linear relationship of bulk salinity and thickness as Cox and Weeks (1974) hypothesized on a limited set of cores, but the slope of the linear relationship in the model is steeper than that proposed by Cox and Weeks (1974).
Our results show that the largest interannual variations of salinity occur at the surface of first-year ice and are caused by rain, surface melt, and flooding.In contrast, the lower 80 % of the salinity profile of first-year ice are similar to each other despite being forced with reanalysis data taken from different locations and years.The multiyear ice profiles vary depending on the ice thickness at the onset of growth and become more similar over the growth season.
We compare the ice thickness, freshwater column, thermal resistance, and total stored energy of the nine 4.5-year simulations of Arctic sea-ice using the three different salinity approaches against each other.Although certain quantities differ by up to 12 % for a specific simulation, on average the differences between the complex salinity approach and the other approaches are below 5 %.The simple approach has a roughly 30 % smaller difference compared to the complex approach than the prescribed approach (Fig. 17) and a roughly 50 % better ratio than the prescribed approach (Fig. 16).
Given that the strong arctic halocline should prohibit strong ice-ocean feedbacks, we expect that fully parametrizing the temporal sea-ice salinity evolution in the Arctic will not have a large effect on sea-ice thermodynamics in climate models.We expect that parametrized-prescribed hybrids, such as that proposed by Vancoppenolle et al. (2009), which parametrizes the evolution of the bulk salinity of the whole ice column and prescribes an empirical salinity profile based on the bulk salinity, will reproduce the dominant thermodynamic effects of the sea-ice salinity evolution.Prescribing age-and thickness-dependent salinity profiles such as those shown in Fig. 14 is also a viable alternative.The multiyear profile of Schwarzacher (1959) underestimates the mean salinity profile because it is based on cores taken from May to September, during which the salinity is lower (Fig. 12).A smooth transition from first-year to multiyear ice can be achieved by linearly transitioning from the first-year to the multiyear profile as discussed in Sect.4.4.Further refinement can be achieved by taking into account the annual cycle (Fig. 12), the ice thickness (Fig. 9), and the sea-ice location.
Comparisons to laboratory and field salinity measurements have shown that the parametrized brine fluxes in SAMSIM are a reasonable approximation of reality.SAM-SIM's semi-adaptive grid is convenient when studying processes which occur close to the ice-atmosphere or ice-ocean boundary, as it avoids numerical diffusion through layer advection in the surface and bottom layers.All dissolved tracers in brine can be easily advected similar to salt, and the gas volume fraction in each layer can be used to compute outgassing and uptake.Thanks to these properties SAMSIM is a valuable tool to study small-scale thermodynamic and other aspects of sea ice which are affected by brine dynamics such as sea-ice biology.

Figure 1 .
Figure 1.Sketch of SAMSIM grid evolution for three top ice layers during snow melt and following surface ablation as explained in Sect.2.3.

Figure 5 .
Figure 5. Brine fluxes of the complex flushing parametrization resulting from meltwater formation at the surface as described in Sect.2.4.2.The horizontal fluxes f h transport heat and salt to the lowest layer directly via cracks in the ice, while the vertical fluxes f v advect heat and salt from layer to layer.ζ is the freeboard of the ice and z melt is the depth of the meltwater.

Figure 6 .
Figure 6.Temperature and bulk salinity evolution of the idealized flushing experiment using the default model set-up (experiment setup in Sect.3, model set-up in Table1).Plot background is grey.The 3 and 7 g kg −1 contour lines are included.

Figure 7 .
Figure 7. Salinity evolution of the idealized melting experiments in which one specific parameter or setting has been changed from the default values (default model results shown in Fig.6, experiment description can be found in Sect.3, default settings are listed in Table1).The 3 and 7 g kg −1 contour lines are included.(a) Gravity drainage is included, which is otherwise disabled in the experiment.(b)Theratio of horizontal to vertical hydraulic resistance β is 0.2 instead of 1.0.In (c) the ratio of horizontal to vertical hydraulic resistance β is 5 instead of 1.0.In (d) the vertical spatial resolution z 0 is 2 mm instead of 1 cm, and in (e) the vertical spatial resolution z 0 is 5 cm instead of 1 cm.

Figure 8 .
Figure 8. Salinity evolution of the (a) complex, (b) simple, and (c) prescribed salinity approach for one of the nine Arctic simulations forced with ERA-interim data from 75 • N, 145 • W. The salinity approaches are described in Sect.2.5, and the model set-up is described in Sect.4.1.The simulation time (x axis) begins on 1 July 2005.The dashed line marks the snow-ice boundary.The water surface is at z = 0.

Figure 9 .
Figure 9.The vertically integrated vertical bulk salinity as a function of ice thickness for all reanalysis forced runs as described in Sect.4.1.Each grey dot represents a 12-hourly snapshot.(a) Contains all 15 years of first-year ice and (b) contains all 21 years of multiyear ice in grey.Of all nine simulations, a single simulation is plotted in black (80 • N, 90 • E) to enable tracking the evolution over time.The blue curve in (a) is the empirical relationship for first-year ice published byKovacs (1997) for ice up to 2 m.The red dashed lines mark the empirical linear relationships found byCox and Weeks (1974) for growing (upper lines) and melting Arctic ice (lower line).

Figure 10 .
Figure10.Time-averaged and vertically normalized salinity profiles from first-year ice cores (described in Sect.4.3 and shown in a) and first-year ice from reanalysis forced simulations using the complex brine dynamic parametrizations (b).Both were averaged from January to March (1-3), April to May (4-5), and over June (6).
Figure11.May to September mean of vertically normalized multiyear salinity profiles of reanalysis forced simulations using the complex brine dynamic parametrizations.Schwarzacher 59 refers to the fitted profile ofSchwarzacher (1959), and Barrow cores refers to the multiyear ice cores taken by the Alaska Ocean Observing System from 1999 to 2011(Eicken et al., 2012).The R crit , α spread shows the SAMSIM profile using the two non-default values of the gravity drainage parameters obtained from the optimization process (see Sect. 2.2).Left line: α = 0.000681, R crit = 3.23.Right line: α = 0.000510, R crit = 7.10.The area between the two simulations is shaded in light grey.

Figure 12 .
Figure 12.Monthly mean of vertically normalized salinity profiles of reanalysis forced simulations using the complex brine dynamic parametrizations as described in Sect.4.1.The simulations were split into annual cycles beginning in September (month 9) and sorted into 15 years of first-year ice (a and b) and 21 years of multiyear ice (c and d).The corresponding ice thickness of the monthly means are shown in Fig. 13.

Figure 13 .
Figure13.The white columns show the thickness of all monthly mean salinity profiles shown in Fig.12.The black columns represent only first-year ice which evolves into multiyear ice the following year.To be included in the monthly average ice must be present, meaning that model output of ice-free water with an ice thickness of zero is excluded from the mean.
along with the mean SAMSIM profiles.The fitted first-year ice profile is S bu, fy (z) = z a + bz + c (17) for a = 1.0964, b = −1.0552,and c = 4.41272, and the fitted multiyear ice profile is S bu, my (z) a = 0.17083, b = 0.92762, and c = 0.024516.

Figure 15 .
Figure 15.Vertically normalized salinity profiles of the reanalysis forced simulations (described in Sect.4.1) using the complex salinity parametrizations on 1 November (a and c) and 1 April (b and d).First-year ice (a and b) and multiyear ice (c and d) are shown separately.The grey lines are the individual model realizations and the black line is the average overall profiles.
Figure17.Time-integrated absolute differences of the ice thickness, freshwater column, thermal resistance R th , and enthalpy H of simulations using the simple and prescribed SAMSIM salinity approach compared to simulations using the complex approach (details in Sect.5).The absolute differences were calculated separately for each of the nine reanalysis forced simulations over 4.5 years.Each dot shows the ratio of a specific simulation, while the lines show the mean overall runs and quantities.

Table 1 .
Default model settings and free parameter values of salinity parametrizations.
Sketch of snow melt by snow-to-slush conversion as described in Sect.2.3.1.Snow-to-slush conversion occurs when the liquid fraction exceeds l, max as shown in the left sketch.A slush layer of thickness B is formed, which is instantaneously added to the top ice layer.A is the thickness lost by snow-to-slush conversion.The top ice layer thickness increases by B while the snow layer thickness is reduced by A + B. The white, blue, and grey areas represent the solid, liquid, and gas volume fractions of the model layers.The combined solid and liquid volumes of the snow and top ice layer are conserved during the conversion.

Table 2 .
Horizontal and vertical fluxes of the thought experiment detailed in Sect.2.4.3 and shown in Figure