Version 1 of a sea ice module for the physics based, detailed, multi-layer SNOWPACK model

multi-layer SNOWPACK model Nander Wever1,2,3, Leonard Rossmann4, Nina Maaß5, Katherine C. Leonard1,2,6, Lars Kaleschke4, Marcel Nicolaus4, and Michael Lehning2,3 1Department of Atmospheric and Oceanic Sciences, University of Colorado Boulder, Boulder, CO, USA. 2CRYOS, School of Architecture, Civil and Environmental Engineering, EPFL, Lausanne, Switzerland. 3WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland. 4Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. 5University of Hamburg, Germany. 6Cooperative Institute for Research in Environmental Science (CIRES), University of Colorado Boulder, Boulder, CO, USA. Correspondence: NANDER WEVER (nander.wever@colorado.edu)

snow. However, for deriving sea ice thickness and snow depth from satellite remote sensing, brightness temperatures of snow are also required. These depend strongly on snow stratigraphy and snow properties. Several studies used the snow thermodynamics model SNTHERM, either using the internal sea ice module or coupled to a sea ice model (Jordan et al., 1999;Nicolaus et al., 2006;Chung et al., 2010;Fuller et al., 2015) to improve the thermodynamical upper boundary condition for sea ice, and assessing the snow microstructure on sea ice.
Rather than improving the representation of snow in sea ice models, this paper presents a new sea ice module developed for the physics based, detailed, multi-layer snow cover model SNOWPACK. Our motivation is that the SNOWPACK model has a long development history with a focus on accurately representing physical processes in the snow cover. This includes a detailed representation of snow microstructure as well as liquid water flow processes (Bartelt and Lehning, 2002;Lehning et al., 2002a, b;Wever et al., 2015Wever et al., , 2016. The model has previously been successfully applied in Polar regions. For example, 10 for the Greenland Ice Sheet (Steger et al., 2017), it was found to provide accurate simulations of water flow and refreezing processes. An application to the Antarctic Plateau (Groot Zwaaftink et al., 2013) showed good agreement for new snow density and temperature profiles in a drifting snow dominated environment.

Methods
The SNOWPACK model has a long development history regarding snow processes. The model calculates the snow energy 15 balance, the resulting temperature distribution, snow settling (densification), liquid water flow and the evolution of snow microstructure. The snow microstructure is described by four parameters: grain size, bond size, sphericity and dendricity. Governing state equations describe the time evolution of these parameters given snow temperature, density, liquid water content (LWC), etc. The full model description is presented in Bartelt and Lehning (2002); Lehning et al. (2002a, b).
The basic model structure of SNOWPACK is congruent with the mushy-layer theory for sea ice (Hunke et al., 2011;Turner 20 and Hunke, 2014). The model is 1-dimensional, with an arbitrary number of vertical layers of arbitrary depth. Typical layer depth, however, is 1-2 cm. Each layers total volume is subdivided into a part consisting of respectively ice, water and air.
Henceforth these layers are called "snow" layers. Note that SNOWPACK also considers soil as a category, which is hereafter ignored in the context of sea ice.
The sea ice extension of the model provides modifications to the model code to include the effect of salinity on thermal 25 properties and liquid water flow. Furthermore, ice growth and melt at the bottom of the domain is assessed. Flooding is considered to occur only in one-dimension, and lateral advection of liquid water is ignored (Maksym and Jeffries, 2000).
Below, we detail the modifications of the SNOWPACK model to make it suitable for sea ice simulations. The modifications are implemented in the main version of SNOWPACK, and sea ice specific settings are only needed in the configuration files for the model. The code base of the SNOWPACK model is the same and future developments in other parts of the code will 30 also be directly available for the sea ice version.

Heat equation and thermodynamics
SNOWPACK solves the heat equation using finite elements, as described in Bartelt and Lehning (2002), allowing to have a simulated snow surface temperature as well as a temperature at the bottom of the snow column. Each snow layer, called element, has two nodes with the adjacent elements. Several modifications were necessary to take into account the effect of salinity on thermodynamical properties. 5 First, it is assumed that all salt is concentrated in the liquid water in the brine pockets and that the volumetric ice content is free of salt, such that one can write: Here, S is the bulk salinity (g kg −1 ), S b is the brine salinity (g kg −1 ) and θ is the volumetric liquid water content (m 3 m −3 ).
The melting point T m (K) of each snow element can then be expressed as a function of brine salinity by the commonly used 10 relationship (Assur, 1960): where µ is a constant (0.054 K (g kg −1 ) −1 ) and T 0 is the melting point of fresh water (273.15 K).
To solve the heat equation, we assume equilibrium between the element temperature T e (K) and the brine melting point T m . When the ice is heating (cooling), brine is supposed to melt (freeze) instantaneously, in order to maintain T e = T m . This 15 is achieved by feeding back the energy associated with the phase change as a source/sink term in the heat equation (see Eq. 3 in Bartelt and Lehning (2002)). Note that the latent heat released (used) when water freezes (melts) increases (decreases) the temperature locally, thus counteracting the temperature change. In turn, the refreezing (melting) ice would impact brine salinity and thereby the melting temperature. These competing processes slow down the convergence in the solver for the heat equation, which can be mitigated by providing an improved estimate of the melting temperature which satisfies: (i) the 20 condition provided by Eq. 1 for the new LWC, (ii) Eq. 2 for the new melting temperature given the new brine salinity, and (iii) the change in ice content for the given deviation of the new element temperature from the melting point, by algebraically solving these three conditions with the three unknowns.
The bottom node of the domain represents the interface between sea ice and ocean and its temperature is prescribed as a Dirichlet boundary condition using the freezing temperature of the ocean water, which is determined by the prescribed ocean 25 salinity. Typically heat is advected in the ocean from below, referred to as the ocean heat flux. We determine the net energy loss or gain at the bottom node, given the prescribed ocean heat flux and the internal heat flux in the lowest sea ice element.
This net energy is translated into bottom ice growth or melt.
An uncertainty for ice growth is the ice porosity of the newly formed ice. We apply a similar approach to the one presented in Griewank and Notz (2013), where an ice content threshold is defined. When the lowest element has an ice content below 30 this value, the net energy is used to create a new ice element with a brine salinity equal to ocean salinity (Vancoppenolle et al., 2010). Otherwise, the bottom element grows in length. We set a threshold of 0.99 m 3 m −3 , which we also prescribe as the maximum allowed ice content of a single layer. Mass loss is applied by reducing the element length.

Brine dynamics
The SNOWPACK model is equipped with two water transport schemes: a simple bucket scheme and a solver for the full Richards equation for transport in porous media (Wever et al., 2014). Below, we will outline how these water transport schemes were modified for sea ice.
A general modification is that the density of liquid water (ρ) is adjusted for salinity as: where ρ w is the density of fresh liquid water (1000 kg m −3 ), β is a salinity coefficient, approximated as 0.824 kg 2 g −1 m −3 (Massel, 2015, Appendix A).

Bucket Scheme for simple water transport
One of the simplest approaches to simulate liquid water flow through porous media is the so-called bucket scheme, which is 10 also available for the sea ice module in SNOWPACK. The scheme is described in Wever et al. (2014). For sea ice, the bucket scheme has been extended with a simple description of flooding. The pore space of the part of the simulated domain which is located below sea level is filled with liquid water with ocean water salinity S o . The salinity distribution of the sea ice is described by a fixed salinity profile.
The position of the sea level inside the model domain z sl (m) is determined from the isostatic balance: where SWE is the snow water equivalent, defined as the sum over all elements N of the mass of each element j in the model: where θ i,j and θ j is the volumetric content (m 3 m −3 ) for layer j of ice and water, respectively. ρ i is the density of ice

Richards equation for explicit water transport
The mixed form of the Richards equation reads: where t is time (s), z is the vertical coordinate (m), κ is the permeability (m 2 ), η is the viscosity (m s −2 ) and s is a source/sink term (m 3 m −3 s −1 ). The pressure p can be considered the sum of water potential and gravity potential: where h is the pressure head (m), g is the gravitational acceleration (m s −1 ), ρ the density of the flowing liquid (kg m −3 ) and γ is the slope angle. Note that SNOWPACK can be used in sloped terrain. For completeness of the model description, we keep this term, although for sea ice γ is obviously 0. The permeability κ is replaced by the hydraulic conductivity K, which relates to κ via: where µ is the dynamic viscosity of water (0.001792 kg (m s) −1 ).
A critical assumption in the typical application of the Richards equation is that it is assumed that both g and ρ are constant in time and z, and consequently can be eliminated from the equation. Due to salinity variations in sea ice, variations of density of the flowing liquid can occur and are actually considered the driving mechanism in the temporal and spatial evolution of the 15 salinity of sea ice.
Therefore, we rewrite the Richards equation for sea ice by considering ρ as a function of z and only eliminating g, arriving at: As outlined in Wever et al. (2014), the Richards equation is implemented in SNOWPACK by the discretization proposed by 20 Celia et al. (1990). The backward Euler approximation in time coupled with a simple Picard iteration, as shown in Eq. 14 of Celia et al. (1990), becomes, for Eq. 9: where n and m denote the time and iteration level, respectively. Here, we have used to chain rule to write: 25 The last term on the right hand side expresses the liquid water flow driven by density differences. After applying a Taylor expansion to Eq. 15 of Celia et al. (1990), and defining δ m ≡ ρh n+1,m+1 − ρh n+1,m , we arrive at: Finally, Eq. 16 in Celia et al. (1990) becomes: After applying the standard finite difference approximation in space, Eq. 16 in Celia et al. (1990) translates into (i denoting the spatial coordinate): The system of equations described by Equation 14 forms a tri-diagonal matrix. As in Wever et al. (2014), the function DGTSV from the LAPACK library (Anderson et al., 1999) is called to compute the solution. When LAPACK is not available, or 15 not selected on compile time, the Thomas algorithm is used as the implemented default alternative, which does not depend on external libraries. However the Thomas algorithm is not the preferred option as in contrast to DGTSV, it lacks partial pivoting and may suffer from numerical instabilities.

Salinity
The governing equation in 1-dimension for concentration describes the change in salinity as a combination of a diffusion and 20 advection process: Where D is the diffusion coefficient (m 2 s −1 ), considered here as a first approximation independent of temperature, taken as 10 −10 m 2 s −1 (Poisson and Papaud, 1983). q denotes the liquid water flux (m s −1 ).
Here, taking f = 1 results in a fully implicit scheme, whereas taking f = 0.5 corresponds to the Crank-Nicolson scheme (Crank and Nicolson, 1996). The equation is solved after every time step for liquid water flow, assuming that the water flux q is constant with time and can be referenced with the time level n. Furthermore, for LWC θ, both θ n , as well as θ n+1 are 10 known from solving Eq. 14. The domain in SNOWPACK is typically nonuniform and the spatial discretizations in Eq. 16 for the nonuniform grid are based on Veldman and Rinzema (1992). The system of equations described by Equation 16 forms a tri-diagonal matrix, similar to and solved in the same way as the equation for liquid water flow (Eq. 14).
SNOWPACK by default uses the Crank-Nicolson scheme. The fully implicit scheme is only first order accurate, whereas the Crank-Nicolson scheme is second order accurate. Both schemes are unconditionally stable and suffer only minimal numerical 15 diffusion for advection. As with many other common schemes, the advection part does not perfectly conserve sharp transitions.
The Crank Nicolson scheme is, in spite of being unconditionally stable, prone to spurious oscillations in the solution. To choose adequate time steps, w apply the CFL condition, typically required for stability of an explicit scheme, also for the Crank-Nicolson scheme: Note that if the CFL condition is violated, the time step is reduced and the last time step for the Richards equation is also repeated with the reduced time step.

Boundary Conditions
The boundary conditions for the Richards equation (Eq. 6) are determined by a Neumann boundary condition at the top, consisting of the top water flux from rain, evaporation or condensation and a Dirichlet boundary condition by prescribing the 25 bottom pressure head. The pressure head at the bottom of the sea ice corresponds to the water pressure at that depth, which equals the sea level (Eq. 4).
The boundary conditions for the advection terms of the advection-diffusion equation (Eq. 15) are prescribed as a Neumann boundary condition, with fresh water at the top, and ocean salinity at the bottom. For the diffusion term, we implemented no-flux boundary conditions. These are derived by considering the central differences scheme for the diffusion term on a nonuniform grid, determined according to: A no-flux boundary condition can be achieved by forcing gradients over the boundaries to be 0, such that either the left (upper boundary) or right (lower boundary) term in the numerator of Eq. 18 vanishes. 5

Hydraulic Properties
For solving the Richards equation and the salinity transport equation, several parameters which depend on the snow and ice microstructure need to be specified. For saturated hydraulic conductivity (K sat ), we define elements with a porosity (i.e., 1−θ i ) larger than 0.25 as snow, and smaller than 0.25 as ice.
For snow elements, a formulation based on Calonne et al. (2012) is typically used (see Wever et al. (2014)): where r es is the equivalent sphere radius (m). Note that Eq. 14 in Wever et al. (2014) erroneously shows a factor 0.75 (which corresponds to r es being a grain diameter) instead of 3.0.
For sea ice, the saturated hydraulic conductivity (K sat ) is based on Golden et al. (2007): 15 In unsaturated conditions, the van Genuchten-Mualem model (Mualem, 1976) is used to relate the hydraulic conductivity in saturated conditions (Eq. 19 and 20) to unsaturated conditions. Furthermore, in unsaturated conditions, the van Genuchten model is used for the water retention curve, which describes the relationship between capillary suction and LWC (van Genuchten, 1980). The coefficients in this parameterization of water retention in snow is based on the work by Yamaguchi et al. (2012) for snow. We extend this parameterization independent of porosity, by absence of any information of water reten-20 tion in sea ice as a function of salinity and LWC. It has a relatively small impact, as the largest part of the sea ice is below sea level and thus saturated.

In-situ Buoys
We apply the sea ice version of SNOWPACK to snow and sea ice properties measured from two pairs of Snow Buoys and one in prep). We construct a time series by averaging the four snow depth sensors of each buoy. In addition to snow depth, the Snow Buoys measure barometric air pressure and air temperature.
Each IMB consists of a 4.8 m long thermistor string, with a vertical sensor spacing of 0.02 m and provides sea ice temperature as well as interfaces from snow, sea ice, and sea water. The instruments are described by Jackson et al. (2013). It turned out to be good practice to co-deploy Snow Buoys and thermistor string IMBs in order to observe snow depth and sea ice properties at the same time. The full data set of all these Lagrangian observations is available from http://www.meereisportal.de (Grosfeld et al., 2015).
The two selected Snow Buoys (Snow Buoys 2014S12 and 2016S31) have long time series and cover different trajectories, as shown in Fig. 1. Unfortunately, IMB 2014T9 collocated with Snow Buoy 2014S12 stopped transmitting data soon after deployment. IMB 2016T41, collocated with Snow Buoy 2016S31 measured for almost the same period as the Snow Buoy. 5 However, comparisons have to focus on the sea ice properties excluding the snow cover on top, because the IMB was deployed directly onto the sea ice surface and the thermistor chain does not measure the snow cover properties, except for the lowest few cm. During the last 18 hours no valid snow depth data were transmitted and the drift speeds were relatively high, suggesting that it is not a transmission or data logger failure, but rather an indication that flow deformation or breakup is the likely cause of the loss of the Snow Buoy. We consider this particular Snow Buoy due to its long time span, even though a collocated IMB data set is not available.

Initial Conditions
To start each simulation, a description of the initial sea ice state is required. Upon installation of each Snow Buoy, the ice thickness, snow thickness and freeboard was determined. For simulations of these Snow Buoys, we distinguish three categories: (i) sea ice below sea level (ice thickness minus freeboard), (ii) sea ice above sea level (freeboard) and (iii) snow.
For the part of the sea ice below sea level, the volumetric ice content θ i was fixed to 0.9, and the volumetric water content 25 θ w was subsequently calculated as: This formulation leaves a small volumetric air content which can be filled when water refreezes and thereby expands. This is currently required for the stability of the numerical schemes in the SNOWPACK model, but in reality refreezing water would increase the pressure in the brine. The element temperature was initialized by the value recorded by the IMB upon installation. 30 As it takes time for the thermistor string to freeze into the ice and adapt to the surrounding ice temperature, this temperature is mostly representative of the ocean water. The brine salinity was set as the salinity for which the melting point corresponds to the measured temperatures. For the part of sea ice above sea level, the volumetric ice content θ i was also fixed to 0.9, but the remaining space was assigned to air content and the bulk salinity was set to 0 g kg −1 . Snow was initialized with a density of 275 kg m −3 , and a grain size of 0.15 mm. The element temperature was derived from the thermistor string measurements. Finally, the depth of each element was set to 0.02 m.

5
Simulations with the SNOWPACK model require air temperature, relative humidity, incoming shortwave radiation, incoming longwave radiation, wind speed and precipitation. Here, we used ERA5 (European Centre for Medium-range Weather Forecasts ECMFW Reanalysis 5) to provide these parameters to drive the simulations. For each timestep and location, the simulated weather at the closest grid point in the ERA5 model was taken.
SNOWPACK has the possibility to either use a precipitation time series as input to determine when snow fall occurs, or to 10 use a time series of snow depth to interpret increases in measured snow depth as snow fall events when simulated snow depth is below measured snow depth (Lehning et al., 1999). In order to use the latter method for the Snow Buoys and base the mass balance on Snow Buoy data, a layer can be marked in the simulations and tracked throughout the snow -sea ice continuum. By marking the layer that corresponds to the reference level for the snow depth measurements, the measured snow depth can be tracked relative to this marked layer. The output routines of the model have been adapted accordingly to reference the output 15 to either the sea level, or to the marked reference layer.
The ocean heat flux determines the ice mass balance at the bottom of the ice. Its value can be highly variable and dependent on ocean conditions below the sea ice (Ackley et al., 2015). From that study, we use a value of 8 W m −2 , unless otherwise noted. Fig. 2 shows an example of model behaviour when an initially dry and fresh ice layer with a thickness of 1.58 m, consisting of 94% ice and 6% air expressed as volumetric content, is positioned in ocean water with a salinity of 35 g kg −1 . The positive pressure head at the bottom of the sea ice corresponds with the pressure exerted by the displaced water. As a consequence, saline water enters the ice matrix (Fig. 2a,c), until it is in equilibrium with the sea level. As the pressure difference of the liquid 25 water inside the ice matrix and the surrounding ocean water is decreasing, the salt influx rate decreases over time (Fig. 2e).

Explicit Brine Dynamics Example
The brine salinity corresponds to the brine salinity of ocean water (35 g kg −1 , see Fig. 2d), corresponding to a bulk salinity of 1.65 g kg −1 (see Fig. 2c). The added mass to the sea ice (Fig. 2b) causes the ice to sink deeper inside the ocean water, decreasing the freeboard (Fig. 2f).  Fig. 3a and 3b show the simulated temperature of the snow-sea ice system for simulations driven by in-situ measured snow depth and ERA5 precipitation, respectively. Dashed lines denote the reference level, i.e., the snow-sea ice interface as determined upon installation of the Snow Buoy. The sea level as calculated from hydrostatic balance is indicated by the solid line.

Temperature validation
The sea level stays below the snow-sea ice interface, i.e., freeboard is positive in the simulations during the whole simulation  Fig. 3c shows the measured sea ice temperatures from the corresponding IMB. Note that for this IMB, the thermistor string does not extend above the initial 2 cm snow layer on top of the sea ice, such that the time evolution of the snow cover is not recorded by the IMB. We find that the IMB confirms the strong cooling of the sea ice during the austral winter months, as found in the simulations, as well as the near-surface warming to the melting point of fresh water shortly before the last transmission by the buoy. This could indicate an overestimation of the surface energy balance by the model, but also by an underestimated new snow density, and thus underestimated thermal conductivity, by the model.
The thermistor string is also used to determine the heat capacity by heating the string for 1 or 2 minutes and looking at the temperature response. By combining the absolute temperatures and the heating rates, the ice-ocean interface has been determined as shown in Fig. 3c. The IMB data confirm the modelling result that the strong negative energy balance at the top of the snow-sea ice system during austral winter has not resulted in an ice thickness increase. The decrease in ice thickness in December 2016 is not reproduced by the model. The trajectory of the buoy shows a marked change in drift direction (Fig. 1), changing from a northward to a southward drifting during this period. This change of direction may have been accompanied by an intrusion of warm ocean water below the sea ice and an increased ocean heat flux. 5 Fig. 3a and 3b show the simulation driven by in-situ snow depth measurements, and by precipitation time series from ERA5, respectively. The in-situ measured snow depth on the Snow Buoys qualitatively shows a good agreement in terms of timing of precipitation events and the general evolution of the snow cover on the sea ice. However, the ERA5 precipitation leads to accumulations in the austral winter time, that are not present in the snow depth time series. This is not necessary a bias in ERA5, because snow erosion by wind can keep the snow depth constant over extended periods of time. Due to this discrepancy, 10 the total snow depth is slightly overestimated by the ERA5 input at the end of the simulation.
LWC, bulk and brine salinity, as well as the bottom salt flux for the simulations driven by measured snow depth are shown in Fig. 4. LWC (Fig. 4a) shows a strong reduction in austral winter due to the freezing brine. In December, occasional surface meltwater from the snow percolates down in the sea ice. Combined with rising temperatures, a high LWC layer forms around and just below sea level. Fig. 4b shows that the bulk salinity of the sea ice hardly changes over the course of the simulation, 15 whereas the brine salinity (Fig. 4c) clearly shows a relationship with the temperature. This reflects the prescribed thermal equilibrium between brine and the ice, assuming that the brine is at melting temperature. Fig. 4d shows that the added weight of snow accumulation pushes the sea ice deeper into the ocean, increasing the pressure head at the bottom of the sea ice. This leads to an influx of saline water.

20
In Fig. 5, we show an example of a simulation with a prescribed, time-invariant salinity profile. The bulk salinity profile is fixed between the bottom and top of the sea ice, with a value of 12 g kg −1 . The vertical distribution of bulk salinity is then prescribed by a sine function with an amplitude of 8 g kg −1 . The snow is assigned a constant salinity of 1 g kg −1 . The water transport is described by the bucket scheme, which is the only scheme that works with prescribed salinity profiles.
The salinity in snow causes the presence of some amount of liquid water in the snowpack, such that the brine salinity satisfies 25 the melting point temperature. We also find that the high bulk salinity prescribed near the snow -ice interface requires high values of LWC to arrive at the brine salinity that would correspond to the ice temperatures. The latent heat required for melting of ice to achieve this, results in lower temperatures of the sea ice. Overall, the simulations compare well to simulations with explicit salinity transport, in terms of ice thickness and temperature distributions.  higher w.r.t. to the snow depth sensor, as indicated by the solid line. This is congruent with a negative freeboard. It indicates that significant flooding occurred, as also evidenced by the LWC in Fig. 6b and associated high bulk salinity (Fig. 6c). Fig. 7b shows that much of the sea water that flooded the snow refroze, adding substantial ice mass. Fig. 6e shows that flooding also leads to a strong influx of salt to the snow-sea ice system. The flooding saturates the snow, which has significant pore space compared to the sea ice. Therefore, snowfall events of similar magnitude have different effects Bulk Salinity (g kg -1 ) (c) Figure 5. Example simulation for buoy 2016S31, where measured snow depth was used to derive the precipitation events, and the bulk salinity profile was prescribed, for temperature (a), LWC (b) and bulk salinity (c). In (b), dry snow is colored grey. The depth on the y-axis is defined relative to the sea ice-snow interface, as determined upon installation (dashed line). The solid line denotes the sea level as determined by the simulations.

Flooding and Superimposed Ice Formation
Note that the flooding, as depicted in simulations with the Richards equation coupled to the transport equation, is governed by the hydraulic conductivity of ice. In cold ice, hydraulic conductivity can be so low that negative freeboard remains for extended period of times, even without flooding. On the other hand, flooding may also be triggered by deformation and cracking of the sea ice, combined with lateral flow effects. However, Maksym and Jeffries (2000) have shown that the simpler assumption (i.e., negative freeboard will trigger flooding) can already yield satisfying results. In the model, the maximum ice contents is fixed 5 to 0.99 m 3 m −3 , such that hydraulic conductivity is typically large enough for instantaneous flooding.

Forced Warming and Cooling
We    accumulating on below-freezing sea ice with low porosity (see Fig. 8b). This process is visible between April 14 and April 27 and leads to superimposed ice formation ( Fig. 9b and 9c). Upon continued warming, fresh water starts flushing the ice, leading to a rapid reduction of bulk and brine salinity (Fig. 8c,d). and no incoming shortwave radiation. The ocean heat flux was set at 8 W m −2 . We find a consistent cooling and increase in brine salinity (Fig. 10d). This is achieved by freezing the liquid water, as shown by the decrease in LWC (Fig. 10b). The bulk salinity remains approximately constant (Fig. 10c).
An increase in brine salinity also increases the density of the brine. This may lead to flushing of the sea ice, when the heavy 10 brine moves downward and is replaced by lighter ocean water. However, in our simulations, there is only a very small outflow of salinity at the bottom (( Fig. 10e). This result shows that without a description of the convective processes in sea ice resulting from cooling (Griewank and Notz, 2013), the salinity depletion found due to cooling is strongly underestimated by this model approach.

15
A final test is run by starting with only 2 cm of ice, with constant atmospheric conditions to simulate thin ice evolution. The atmospheric conditions were set as -10 • C air temperature, 100% relative humidity, no wind speed, no incoming solar radiation and a constant incoming longwave radiation of 230 W m −2 . The ocean heat flux was set to 0 W m −2 .
The simulations were run for 1 month, which resulted in approximately 50 cm ice growth (Fig. 11). The relatively warm sea ice compared to the air temperature results in a latent heat flux directed to the atmosphere, even though relative humidity 20 is 100%. The evaporation at the top of the sea ice leads to an outflow of fresh water at the top of the snowpack, resulting in an accumulation of salt near the surface (Kaleschke et al., 2004;Domine et al., 2005, e.g.). The transport of salt from below decreases with increasing ice thickness, as capillary forces are not strong enough anymore to bridge the freeboard.

Conclusions
We introduced a series of modifications to the physics-based, multi-layer SNOWPACK model to allow it to simulate sea ice. 25 This involved modifying the model thermodynamics to account for the varying melting point of ice based on salinity and adding domain restructuring to allow basal ice growth. Water transport through the snow -sea ice system can optionally be described the process is largely underestimated compared to what is known from literature, as convective brine transport is, thus far, not described by the model.
The modifications allow the snow microstructure descriptions developed in the SNOWPACK model to be applied for sea ice conditions as well. The model is able to simulate the temporal evolution of snow density, grain size and shape and snow wetness over the life span of an ice floe. We find abundant depth hoar layers and melt layers, as well as superimposed ice 5 formation due to flooding and percolation.
Driving the simulations, particularly estimating the ocean heat flux from below, is difficult due to limited forcing data from polar regions. Simulations based on Snow and IMB Buoy data demonstrate the importance of such remote data collection systems for modelling. The detailed snow microstructure evolution can be used for remote sensing retrieval algorithms to assess snow depth and ice thickness from space. For this, a coupling to the Snow Microwave Radiative Transfer model SMRT 10 (Picard et al., 2018) is foreseen.
Code and data availability. The SNOWPACK model and the MeteoIO meteorological preprocessing library (Bavay and Egger, 2014) needed to run SNOWPACK are available under a LGPLv3 license under https://models.slf.ch. The source code of the version used for the simulations presented in this study corresponds to revision 2380 of MeteoIO (https://models.slf.ch/svn/meteoio/trunk) and revision 1751 of SNOWPACK (https://models.slf.ch/svn/snowpack/branches/dev). The input and configuration files for the example simulations in this study are available as an Online Supplement. The website https://niviz.org/ can be used to visualize the SNOWPACK output files.
Author contributions. NW developed the model code and performed the simulations, LR performed code testing, LR, NM, KL, LK, MN and ML contributed to model conceptualization and design. NW prepared the manuscript with contributions from all co-authors.