# Snow topography on undeformed Arctic sea ice captured by an idealized “snow
dune” model

Predrag Popović  All correspondence should be directed to Predrag Popović
<EMAIL_ADDRESS>The University of Chicago, The Department of the Geophysical
Sciences
Justin Finkel  The University of Chicago, Committee on Computational and
Applied Mathematics
Mary C. Silber  The University of Chicago, Department of Statistics and
Committee on Computational and Applied Mathematics
Dorian S. Abbot  The University of Chicago, The Department of the Geophysical
Sciences

###### Abstract

Our ability to predict the future of Arctic sea ice is limited by ice’s
sensitivity to detailed surface conditions such as the distribution of snow
and melt ponds. Snow on top of the ice decreases ice’s thermal conductivity,
increases its reflectivity (albedo), and provides a source of meltwater for
melt ponds during summer that decrease the ice’s albedo. In this paper,
Resub.19we develop a simple model of pre-melt snow topography that accurately
describes snow cover of flat, undeformed Arctic sea ice on several study sites
for which data was available. The model considers a surface that is a sum of
randomly sized and placed “snow dunes” represented as Gaussian mounds. This
model generalizes the “void model” of Popović  _et al._ (2018) and, as such,
accurately describes the statistics of melt pond geometry. We test this model
against detailed LiDAR measurements of the pre-melt snow topography. We show
that the model snow-depth distribution is statistically indistinguishable from
the measurements on flat ice, while small disagreement exists if the ice is
deformed. We then use this model to determine analytic expressions for the
conductive heat flux through the ice and for melt pond coverage evolution
during an early stage of pond formation. We also formulate a criterion for ice
to remain pond-free throughout the summer. Results from our model could be
directly included in large-scale models, thereby improving our understanding
of energy balance on sea ice and allowing for more reliable predictions of
Arctic sea ice in a future climate.

##  I Introduction

As a hallmark of climate change and a major component of the Arctic
environment, Arctic sea ice retreat has garnered much scientific and media
attention (Perovich and Richter-Menge, 2009; Stroeve  _et al._ , 2007).
Predicting the rate of sea ice loss is a challenge. The difficulty lies in the
fact that sea ice evolution is controlled by many processes that operate on
scales ranging from sub-millimeter to tens of kilometers (Holland and Curry,
1999). Possibly the most notable among these processes is the interaction of
ice with fluxes of energy coming from the environment mediated by detailed
conditions on the ice surface. The presence of snow, water, dirt, Resub.41or
black carbon on the ice surface can drastically change the rate of absorption
of solar radiation or the rate of thermal conduction during winter growth
(Perovich, 1996). Finding ways to reduce the complexity of modeling ice
surface conditions is of great importance for accurately determining the ice
energy balance in large-scale models and, therefore, for improving our
understanding of the future of sea ice in a changing climate.

Snow insulates the ice (Yen, 1981; Sturm  _et al._ , 2002), reflects sunlight
(Perovich, 1996), and provides a source of fresh water that collects into melt
ponds (Polashenski  _et al._ , 2017). The net effect of each of these
processes depends on the spatial distribution of the snow cover.
ResubResub.67This is important because snow often exhibits significant spatial
variability that arises from a multitude of bedforms, most notable of which
are meter-scale snow dunes (Mather, 1962; Filhol and Sturm, 2015). Ice covered
with patchy snow will, on average, grow faster than if the same amount of snow
were spread uniformly, since uniform cover insulates the ice more thoroughly
(Sturm  _et al._ , 2002). Second, uniform snow cover also protects the
underlying ice more from solar radiation than a patchy layer, since even a
thin layer of snow can increase albedo significantly (Perovich, 1996).
Finally, as snow starts to melt, melt ponds first appear in the regions where
there is the least snow (Petrich  _et al._ , 2012). For this reason, patchy
snow cover leads to melt ponds covering the ice surface sooner, melting more
ice. In summary, patchy snow cover will lead to both more ice growth during
winter as well as to more ice melt during summer.

Liston  _et al._ (2018) modeled snow on flat ice using a simple statistical
approach. They represented the snow surface by smoothing and rescaling an
initially random height field to match the measured mean snow depth, its
variance, and the horizontal correlation length. Their approach is likely
sufficient for many practical applications where snow topography needs to be
included. However, it produces some essentially unphysical predictions, such
as negative snow depths in cases when snow depth standard deviation is
comparable to the mean. Here, we adopt a somewhat similar approach that is
physically consistent and reproduces the measurements with high accuracy.

As we mentioned above, snow often exists in meter-scale dunes which are the
most prominent features of the snow cover (Mather, 1962; Filhol and Sturm,
2015). In Popović  _et al._ (2018), we considered a simple geometric “void
model” which represented these snow dunes as circles, and compared it to
aerial photographs of melt ponds. These circles had varying size and were
placed on a surface randomly and independently of each other, while melt ponds
were represented as voids between these circles. Our model was able to
reproduce the geometric statistics of melt ponds, such as their size
distribution and the fractal dimension as a function of pond size, over the
entire observational range of more than 6 orders of magnitude in pond area.

In this paper, we generalize the two-dimensional discontinuous “void model” to
a continuous synthetic “snow dune” topography that has a vertical component.
We represent snow dunes as mounds of Gaussian form that have a randomly chosen
horizontal scale, a height proportional to that scale, and that are placed
randomly on the surface. The surface topography is the sum of many such
mounds. Like the “void model,” this topography accurately describes the
horizontal melt pond features, which we take as indirect evidence that it
accurately describes the horizontal snow features. We support this by showing
that our model fits both LiDAR scans of the pre-melt snow topography and melt
pond data using the same typical horizontal mound scale. The fact that the
horizontal scale is so similar is somewhat surprising since the two datasets
were recoded in different years, different times of year, and different ice
types.

The novel prediction of our continuous topography is a vertical snow-depth
distribution. By comparing moments of our model distribution to LiDAR
measurements of snow-depth, we show that our model height distribution is
indistinguishable from the measurements for two occasions when the underlying
ice was very flat. The one LiDAR measurement on slightly deformed ice showed
subtle deviation from our model. Thus, we hypothesize that our model is a
highly accurate representation of snow on undeformed ice, while it becomes
increasingly inaccurate on deformed ice. The close agreement between snow
topography and our model is the main result of this paper. Since our model has
only three parameters that are uniquely determined by the mean, variance, and
the correlation length of snow depth, it appears that it is only necessary to
measure these three quantities to completely characterize the statistics of
snow cover on flat, undeformed ice.

After showing the agreement between our model and snow measurements, we
consider some applications of our model. First, we consider how spatial
variability in snow affects the heat conduction through the ice during winter
ice growth. Sturm  _et al._ (2002) conducted measurements of snow
conductivity, $k_{s}$, on deformed multi-year ice, and found that $k_{s}$
inferred from ice growth is $\sim 2.4$ times higher than $k_{s}$ inferred from
direct measurements. They found that about 40% of the increase can be
explained by snow and ice geometry and that horizontal heat transport,
typically neglected in large-scale models, likely contributes significantly.
To understand whether these conclusions also hold for undeformed ice, we solve
a 3-dimensional heat equation within the ice assuming that the snow cover is
well-described by our “snow dune” model. We develop a simple analytic equation
to determine the heat flux conducted through the ice given a set of physical
parameters. ResubResub.110Contrary to Sturm  _et al._ (2002), our model shows
a negligible horizontal heat transport. Moreover, our model shows a smaller
increase in heat conduction than Sturm  _et al._ (2002). This disagreement
could be due to the fact that our model assumes undeformed ice while the field
measurements of Sturm  _et al._ (2002) were made on deformed ice.

Next, we consider the effect of snow topography on early melt pond
development. Using the fact that our model height distribution is well-fit
with a gamma distribution, we write an analytic evolution equation for melt
pond coverage during an early stage when ice is impermeable. This equation for
melt pond evolution enables us to understand how early-stage melt pond growth
depends on measurable environmental parameters. Moreover, it allows us to
derive a simple condition that distinguishes whether ice will remain pond-free
throughout the summer based on snow depth, variance, density and melt rate.

This paper is organized so that the new results are presented in the main
text, while confirming previous results about melt ponds on our topography and
mathematical details are left for the Supplementary Information (SI). In
section II, we describe our model of the “snow dune” topography. Next, in
section III we describe the analytical properties of our model such as the
dependence of moments and the correlation function on model parameters. Then,
in section IV, we compare our model with the measured snow topography. In
section V, we apply our “snow dune” model to determine the conductive heat
flux through the ice. Next, in section VI, we use our model to investigate
melt pond evolution on flat impermeable ice. Finally, in section VII, we
discuss the implications of our results for large-scale studies and in section
and VIII, we conclude. In SI section S1, we confirm that the synthetic “snow
dune” topography predicts ponds that accurately reproduce the geometric
statistics of real ponds during late summer. ResubResub.132In the subsequent
SI sections, we prove the mathematical results and show that our model
predictions are robust against details such as the shape of the mounds or
mound radius.

##  II The synthetic “snow dune” topography

Figure 1: a) A realization of the “void model.” Figure is taken from Popović
_et al._ (2018). b) Synthetic “snow dune” topography. Red colors indicate
topographic highs, while blue colors indicate topographic lows. The upper
bound on the scale bar, here set to 1, is arbitrary. c) Ponds on a “snow dune”
topography. d) A binarized image of a real melt pond photograph taken on
August 14th during the HOTRAX mission.

In this section we will describe our synthetic “snow dune” topography. The
entire rest of the paper hinges on this topography.

In Popović  _et al._ (2018), we developed a simple geometric model where snow
dunes were represented by overlapping circles of varying size placed randomly
and independently of each other on the ice surface, while melt ponds were
represented as voids between these circles. We then compared the statistical
properties associated with this simple geometric model with observations of
late-summer melt ponds. We found a remarkable agreement between various
statistics of real and model melt ponds: by tuning only two model parameters,
the model matched the measured pond area distribution, the pond fractal
dimension as a function of pond area, and two correlation functions that
describe pond size and connectedness, over the entire observational range of
nearly 7 orders of magnitude in pond area. A realization of this model is
shown in Fig. 1a.

Resub.145Motivated by the success of the “void model” in capturing the
conditions on the ice surface, we set about to generalize it to a
3-dimensional “snow dune” topography that has an additional vertical
component, and can, thus, be compared with a pond-free ice surface. We will
show that, in addition to reproducing melt pond features, the “snow dune”
model also matches the observed snow-depth distribution and correlation
function on flat, undeformed ice.

To generalize the “void model” to a continuous topography, we replace circles
with mounds that have a vertical profile. The mounds have a Gaussian shape

$h(\mathbf{x})=h_{m}e^{-(\mathbf{x}-\mathbf{x}_{0})^{2}/2r^{2}}\text{ , }$ (1)

where $h(\mathbf{x})$ is the height of the mound at a point $\mathbf{x}$ on
the surface, $\mathbf{x}_{0}$ is the location of the center of the mound, $r$
is the horizontal scale, and $h_{m}$ is the maximum height of the mound.
Mounds are placed randomly, i.e., $\mathbf{x}_{0}$ can be anywhere on the
surface with equal probability. The horizontal scale, $r$, is randomly chosen
from an exponential probability distribution, $f_{r}$

$f_{r}(r)=\frac{1}{r_{0}}e^{-r/r_{0}}\text{ , }$ (2)

where $r_{0}$ is the typical scale of the mounds. Resub.156BThis distribution
is the same as the distribution of circle radii in the “void model” of Popović
_et al._ (2018), and we chose it mainly due to its simple form and the fact
that it depends on a single parameter, $r_{0}$. Otherwise, the choice of this
distribution is arbitrary and does not affect the main conclusions of our
model so long as it is not long-tailed. To prevent having unrealistically
narrow and high mounds, we prescribed the height of each mound, $h_{m}$, to be
proportional to its horizontal scale, $h_{m}=h_{m,0}r/r_{0}$, where $h_{m,0}$
is the typical mound height. The topography is then a sum of $N$ such mounds
placed on an initially flat surface

$h_{\text{SD}}(\mathbf{x})=\sum_{i=1}^{N}h_{m,0}\frac{r_{i}}{r_{0}}e^{-(\mathbf{x}-\mathbf{x}_{0_{i}})^{2}/2r_{i}^{2}}\text{
, }$ (3)

where $h_{\text{SD}}(\mathbf{x})$ is the height of the “snow dune” topography
at location $\mathbf{x}$. A realization of this topography is shown in Fig.
1b. Ponds that form on this surface after cutting it with a horizontal plane
are shown in Fig. 1c. This is shown alongside a realization of the “void
model” (Fig. 1a) and a binarized image of real melt ponds (Fig. 1d).

There are three parameters in this model. These are the typical horizontal
mound scale, $r_{0}$, the density of mounds, $\rho$, placed within the domain
of size $L$, $\rho\equiv N\frac{r_{0}^{2}}{L^{2}}$, and the typical mound
height, $h_{m,0}$. The first two of these parameters, $r_{0}$ and $\rho$, also
enter the “void model.”

Resub.155Resub.M1.1Although our model mounds bare some resemblance with real
snow dunes, our model is not accurate on the scale of individual dunes. In
particular, snow dunes are certainly not Gaussian-shaped, they show strong
anisotropy, Resub.156the distribution of dune sizes is likely not exponential,
and dunes sometimes form semi-regular lattices, which violates the random
placement assumption (Kochanski  _et al._ , 2019a; Filhol and Sturm, 2015;
Sharma  _et al._ , 2019; Kochanski  _et al._ , 2019b). In SI subsection S2.3,
we show that aggregate, radially averaged statistical properties of the snow
surface that we consider here do not depend on most of these details. We also
note that the random placement assumption is important for the conclusions we
make, so our model can be useful for large scale studies only if snow surfaces
that strongly violate this assumption do not occur often in nature.

##  III Statistics of the synthetic “snow dune” topography

Here, we describe the statistics of the “snow dune” model that we will use
later. In particular, we first show how the mean, variance, and correlation
length depend on model parameters. We then derive all of the moments of the
topography. Finally, we compare the height distribution of the “snow dune”
topography with a gamma distribution. We derive all of the results in this
section in SI section S2. All of these results follow directly from the
definition of the “snow dune” model we described in section II.

First, we find that the mean, $\langle h_{\text{SD}}\rangle$, and variance,
$\sigma^{2}(h_{\text{SD}})=\langle(h_{\text{SD}})^{2}\rangle-\langle
h_{\text{SD}}\rangle^{2}$, of the topography depend on the typical mound
height, $h_{m,0}$, and the density of mounds, $\rho$, as

$\langle h_{\text{SD}}\rangle=12\pi h_{m,0}\rho\quad\text{ ,
}\quad\sigma^{2}(h_{\text{SD}})=24\pi h_{m,0}^{2}\rho\text{ . }$ (4)

In fact, we can explicitly find every moment of the topography,
$\big{\langle}(h_{\text{SD}})^{n}\big{\rangle}$, by knowing the previous
moments and using the recursion formula (see SI section S2 for derivation)

$\big{\langle}(h_{\text{SD}})^{n}\big{\rangle}=2\pi\rho\sum_{j=1}^{n}{n-1\choose
n-j}\frac{(j+2)!}{j}h_{m,0}^{j}\big{\langle}(h_{\text{SD}})^{n-j}\big{\rangle}\text{
. }$ (5)

Equations 4 also follow from this formula.

Next, we define the height correlation function, $C_{h}(\mathbf{l})$, as

$C_{h}(\mathbf{l})\equiv\frac{\langle
h_{\text{SD}}(\mathbf{x})h_{\text{SD}}(\mathbf{x}+\mathbf{l})\rangle-\langle
h_{\text{SD}}\rangle^{2}}{\sigma^{2}(h_{\text{SD}})}\text{ , }$ (6)

where $\mathbf{l}$ is a horizontal displacement of the topography.
ResubResub.209Below, we will use the boldface notation, $\mathbf{l}$, to
denote the vector displacement and we will use the italic notation, $l$, to
denote the magnitude of this displacement. The correlation function quantifies
how much the surface height at $\mathbf{x}$ is correlated with the surface
height at $\mathbf{x}+\mathbf{l}$. If the surface is isotropic, $C_{h}$ only
depends on the magnitude of the displacement vector, $l$. We can then define
the correlation length, $l_{0}$, as the distance at which $C_{h}$ falls by a
factor of $e$. Beyond $l_{0}$, features on the topography can be considered
approximately uncorrelated. Using these definitions, we find that $C_{h}$ and
$l_{0}$ depend on the horizontal mound scale, $r_{0}$, of an isotropic “snow
dune” model as (see SI section S2 for derivation)

$\displaystyle C_{h}(l)$
$\displaystyle=\frac{1}{24}\int_{0}^{\infty}z^{4}e^{-z-l^{2}/(2r_{0}z)^{2}}\text{d}z\text{
, }$ (7) $\displaystyle l_{0}$ $\displaystyle=r_{0}\xi_{0}\text{ , }$ (8)

where $\xi_{0}\approx 9.3689\text{...}$ is a number that can be estimated to
arbitrary precision by inverting $C_{h}(l_{0})=1/e$ using Eq. 7. The mean,
variance, and the correlation length are all measurable properties of real
topographies. ResubResub.216This means that computing these quantities from
measurements lets us unambiguously choose the model parameters limited only by
the accuracy of those measurements.

Figure 2: Dots represent the height distribution of the “snow dune” topography
for different densities of mounds, $\rho$. Solid lines represent fits using a
gamma distribution. Units of snow depth here are arbitrary. Note that the
LiDAR measurements we consider in section IV range from $\rho\approx 0.2$ to
$\rho\approx 0.5$, a much smaller range than we investigated here.

Since we can find all of the moments of the “snow dune” topography using Eq.
5, we can fully determine its height distribution. However, it is impractical
to do this, and, for practical purposes, we now show that the height
distribution of the “snow dune” topography is well-fit with a gamma
distribution. The two distributions are not the same, as can be seen by
comparing their moments (see SI figure S2). However, they are qualitatively
very similar, and the differences between them arise from arbitrary choices in
our model such as the exponential distribution of radii of the snow dunes, Eq.
2, and the linear relationship between the mound radius and height. We
describe the relationship between the height distribution of the “snow dune”
topography and the gamma distribution in detail in SI section S2.

A gamma distribution has the form

$f_{\Gamma}(h)=\frac{1}{\Gamma(k)h_{0}^{k}}h^{k-1}e^{-h/h_{0}}\text{ , }$ (9)

where $h_{0}$ is a scale parameter, $k$ is a shape parameter, and $\Gamma(x)$
is a gamma function. The two parameters of the gamma distribution, $h_{0}$ and
$k$, are related in a simple way to the mean and variance of the distribution,
$\langle h\rangle$ and $\sigma^{2}(h)$

$k=\frac{\langle h\rangle^{2}}{\sigma^{2}(h)}\quad\text{ , }\quad
h_{0}=\frac{\sigma^{2}(h)}{\langle h\rangle}\text{ . }$ (10)

Both of these parameters are therefore uniquely determined by the typical
height and density of mounds through Eqs. 4. In particular, we find that the
scale parameter, $h_{0}$, scales linearly with the typical height of mounds,
while the shape parameter, $k$, scales linearly with the density of mounds,

$h_{0}=2h_{m,0}\quad\text{ , }\quad k=6\pi\rho\text{ . }$ (11)

In Fig. 2, we show height distributions for synthetic topographies with
varying spatial densities of mounds, $\rho$, with the variance,
$\sigma^{2}(h_{\text{SD}})$, kept fixed. We can see that varying $\rho$
changes the shape of the distribution and that a gamma distribution can fit
the height distribution well for all choices of $\rho$. To quantify the
quality of the fit, we use the Kolmogorov-Smirnov (KS) statistic, i.e., the
maximum distance between the empirical and theoretical cumulative
distributions. The fit is typically considered good if this statistic is below
0.05. We find that for all choices of $\rho$, the maximum distance between the
cumulative “snow dune” height distribution and the best-fit gamma distribution
is below 0.05, and ranges between 0.03 for $\rho=0.01$ to 0.006 for $\rho=1$,
generally decreasing for increasing $\rho$. Resub.223Our “snow dune” model
height distribution compares similarly well with a gamma distribution
according to other metrics such as the Kuiper or the Cramér-von Mises metrics.

Resub.M1In SI subsection S2.3, we compare different variants of our model that
include anisotropy, different dune shapes, and different distributions of
mound sizes. There, we show that moments of the height distribution as well as
the radially averaged correlation function are insensitive to these details.
Resub.164BNamely, these alterations were only able to change the numerical
pre-factors in equations that relate model parameters to measurable
quantities, such as $12\pi$ and $24\pi$ in Eqs. 4 or $\xi_{0}$ in Eq. 8. Since
our conclusions regarding heat conduction and melt pond evolution that we
discuss later in the paper only rely on these invariant statistics, they
should also be robust to these. We caution, however, that certain phenomena,
such as the wind stress, may explicitly depend on specifics such as snow
anisotropy. In such a case, additional parameters in the model describing
these details would be required.

##  IV Measured snow topography

Figure 3: a) LiDAR measurement of the pre-melt snow topography by Polashenski
_et al._ (2012) in early June of 2010. b) A realization of synthetic “snow
dune” topography. The density of mounds was chosen such that the height
distribution of the synthetic topography corresponds well to the measured
snow-depth distribution. c) Dots represent snow-depth distribution estimated
using LiDAR measurements for two different years. Solid lines are gamma
distributions with parameters estimated based on the mean and variance of the
measurements. Dashed lines are the corresponding “snow dune” topography height
distributions, and are nearly indistinguishable from the solid lines. d)
Moments of the height distribution. The horizontal axis is the moment order,
$n$, while the vertical axis is the $n$-th root of the $n$-th moment. Dots
represent moments of the measured snow topography. The thick dashed lines
represent the mean moments across an ensemble of 50 randomly-generated “snow
dune” topographies on the domain of the same size as the measurements. The
colored shadings represent one standard deviation from the mean on that
ensemble.

The synthetic “snow dune” topography is meant to represent the height of snow
relative to the flat underlying ice. Therefore, the height distribution of the
“snow dune” topography should correspond to the pre-melt snow-depth
distribution. Here, we will compare the statistics of the synthetic “snow
dune” topography to detailed LiDAR measurements of pre-melt snow topography
made by Polashenski  _et al._ (2012).

During their field expedition, Polashenski  _et al._ (2012) performed detailed
LiDAR scans of the surface topography within a 100m$\times$200m region on
multiple dates in 2009 and 2010 Resub.246near Barrow, Alaska. During 2009,
Polashenski  _et al._ (2012) monitored two locations separated by about 1 km
from each other, one in the north (2009N) and one in the south (2009S). In
2010 they monitored only one location. In Fig. 3a we show an example of such a
measurement, which shows the height of snow before the start of the melt
season, and compare it to a randomly generated “snow dune” topography (Fig.
3b). Assuming that the underlying ice is flat and the pre-melt ice is fully
covered with snow, the LiDAR-estimated height in Fig. 3a represents the snow
depth plus some reference height. Therefore, to compare with the “snow dune”
topography, we estimated the snow depth from these LiDAR scans by subtracting
the minimum LiDAR elevation from the rest of the scan. We show the snow-depth
distribution for the three LiDAR scans in Fig. 3c and we show the moments of
these distributions in Fig. 3d. We estimate the height correlation function
for these three measurements in Fig. 4. The data from Polashenski  _et al._
(2012) are freely available at http://chrispolashenski.com/data.php.

To see whether the snow-depth distribution of these measurements conforms to
the predictions of our model and a gamma distribution, we inferred parameters
$h_{m,0}$ and $\rho$ using Eqs. 4 and parameters of the gamma distribution,
$k$ and $h_{0}$, using Eqs. 10 based only on measured mean snow depth and snow
depth variance. We show the statistics of measured topographies along with
model and gamma distribution parameters in Table 1. The solid lines in Fig. 3c
show a gamma distribution with parameters chosen in this way, while the dashed
lines show a height distribution of the synthetic “snow dune” topography. We
can see that in all cases, the measured snow-depth distribution agrees well
with our model and a gamma distribution, with a Kolmogorov-Smirnov (KS)
statistic remaining at or below 0.03 in all cases.

Table 1: Statistics of the LiDAR measurements and model parameters inferred from them. The KS statistic was calculated with respect to the gamma distribution. | $\langle h\rangle$ | $\sigma(h)$ | $l_{0}$ | $h_{m,0}$ | $\rho$ | $r_{0}$ | $h_{0}$ | $k$ | KS
---|---|---|---|---|---|---|---|---|---
2009N | 15.2 cm | 7.8 cm | 5.5 m | 2.0 cm | 0.20 | 0.59 m | 4.0 cm | 3.8 | 0.03
2009S | 13.4 cm | 5.4 cm | 5.2 m | 1.1 cm | 0.33 | 0.56 m | 2.2 cm | 6.2 | 0.01
2010 | 13.4 cm | 4.3 cm | 5.8 m | 0.7 cm | 0.53 | 0.61 m | 1.4 cm | 9.9 | 0.03

We can gain a better understanding of both the measured and the synthetic
snow-depth distributions by looking at their moments. Higher order moments
describe the more and more detailed structure of the probability distribution.
This is precisely why we can use these higher order moments to distinguish
between similar, but subtly different distributions. In Fig. 3d, we compare
moments of the measured snow topography and randomly-generated “snow dune”
topographies. Higher order moments depend on the domain size, resolution, and
the particular realization of the topography. For this reason, to get the
moments of the simulated “snow dune” topographies, we created an ensemble of
50 randomly-generated topographies with the same resolution and domain size as
the measurements. For each realization, we then found the root-moments,
$\langle h^{n}\rangle^{1/n}$. Dashed lines in Fig. 3d show the mean and one
standard deviation of the root-moments over this 50-member ensemble. We used
the model parameters shown in Table 1.

For the 2010 and 2009S measurements, all of the measured moments fall squarely
within the range of values obtained with an ensemble of simulated “snow dune”
topographies. In fact, we tested the first 150 moments for these two
measurements, and found that this agreement holds throughout. This means that
realizations of the “snow dune” topography that have exactly the same height
distribution as the measurements are common. We can therefore conclude that
the snow-depth distribution for 2010 and 2009S measurements is
indistinguishable from the “snow dune” height distribution! This agreement is,
however, not observed for 2009N measurements - moments of measured snow-depth
distribution are consistently higher than those of the simulations and no
randomly generated “snow dune” topography has high order moments that match
the measurements. This difference is subtle, which is why we could not observe
it using the KS statistic, and for most practical uses that require only a
snow-depth distribution it is unlikely to be important. Nevertheless, it
constitutes a real observable difference between our model and the 2009N
measurements, and shows that there exist features in real data that our model
cannot capture. We return to explaining this difference below.

Figure 4: a) Height correlation functions. Dots represent measurements, and
the solid lines represent predictions for the “snow dune” model using Eq. 7.
Correlation functions for 2010 and 2009S measurements (yellow and blue dots)
were calculated from raw data. Correlation function for the 2009N measurement
was corrected by removing a long-decay exponential. The raw 2009N data are
shown in panel b. b) Correction for the 2009N correlation function. Dots
represent $C_{h}$ calculated with raw data on a semi-log plot. The two black
dashed lines are exponential fits to the short-decay and long-decay portions
of $C_{h}$. Inset shows $C_{h}$ before (red dots) and after (red dashed line)
removing the long-decay exponential on a linear plot.

Next, we describe the height correlation function, $C_{h}(l)$, for
measurements and the model. We find the LiDAR $C_{h}(\mathbf{l})$ using Eq. 6.
A22Due to prevailing winds, snow on sea ice often shows a strong preferential
orientation, with snow dunes elongated along the direction of the wind
(Petrich  _et al._ , 2012). We observe this in the measured
$C_{h}(\mathbf{l})$, which we find to depend on the direction of the
displacement vector, $\mathbf{l}$. Such anisotropy can be included in our
model by elongating the mounds along a particular axis. However, to keep our
discussion as simple as possible, we chose to keep our model isotropic. Thus,
to obtain a measured correlation function that depends only on the magnitude
of $\mathbf{l}$, we average the measured $C_{h}(\mathbf{l})$ over all
directions of $\mathbf{l}$. Then, we find the correlation length, $l_{0}$, the
distance at which the measured $C_{h}(l)$ falls off by a factor of $e$, and
relate it to the model parameter $r_{0}$ according to Eq. 8 (see Table 1). We
show the measured and the theoretical $C_{h}(l)$ calculated according to Eq. 7
in Fig. 4.

Again, we find that 2010 and 2009S measurements agree very well with model
predictions (Fig. 4a). However, 2009N shows some disagreement with the model.
The height correlation function calculated from raw 2009N data is shown in
Fig. 4b. It first decays quickly, and then the decay slows down, with features
on the surface remaining correlated across the domain. By plotting $C_{h}$ on
a semi-log plot, we find that it is approximately a sum of two exponentials -
one with a short decay-length and another with a long decay-length comparable
to the domain size. We get good agreement with the model after removing this
long decay exponential. To do this, we fit an exponential function to $C_{h}$
for $l>20\text{m}$. We then subtract this exponential from the whole $C_{h}$
and rescale so that $C_{h}(0)=1$. The result is a corrected correlation
function, shown as the red dots in Fig. 4a and the dashed line in the inset of
4b. The result of this correction does not depend strongly on the choice of
the cutoff length for the exponential fit so long as this length is large
enough. The correlation length, $l_{0}$, and the model parameter $r_{0}$
reported in Table 1 correspond to this corrected correlation function.

In Popović  _et al._ (2018) and SI section S1, we calculated similar
correlation functions for melt pond images taken during 1998 and 2005. There,
we similarly observed correlation functions with two decay lengths - one
comparable to the size of melt ponds, and another comparable to the size of
the images. In Popović  _et al._ (2018), we attributed the long-decay length
to differences in ice properties between different ice floes and large-scale
deformation features such as ridges. A match between the model and the data
could only be achieved after removing this long-decay length by a procedure
analogous to the one we described above.

Here, we again suggest that the long-decay length comes from variability in
the underlying properties of the ice. In fact, we believe that all of the
differences between 2009N measurements and the model can be explained by this
variability in the underlying ice properties. Polashenski  _et al._ (2012)
note that the 2009N site had ice that finger-rafted early in the growth season
and contained some lightly rubbled ice, while the 2009S site had very flat
ice. They did not comment on ice conditions during 2010 measurements.
Resub.338To test whether our model is consistent with variable ice properties
across the domain, we created a “snow dune” topography where one half of the
domain had higher $h_{m,0}$ and $\rho$ than the other half, but together both
parts had the same mean and variance as the measured 2009N topography. In this
experiment, the model correlation function obtained a long decay length
comparable to the domain size. Moreover, we found that we can choose the
parameters in the first half of the domain such that all of the moments of the
model height distribution match the moments of the 2009N data. Therefore, our
model is consistent with the 2009N data if we assume that 2009N data contains
regions with different underlying ice properties. We note that even though
there exists an observable difference between the snow-depth distribution of
2009N and our model, it is very small and would likely not significantly
change the ice evolution in a large scale model. Nevertheless, it hints that
highly deformed ice may in fact significantly deviate from our model.

In SI section S1, we compare melt ponds on the synthetic “snow dune”
topography with melt pond data derived from images taken during 1998 and 2005.
These images cover a much larger area than the LiDAR scans we considered here.
Each melt pond image covers an area of roughly 1 $\text{km}^{2}$ as opposed to
$0.02\text{ km}^{2}$ covered by the LiDAR scans. As in Popović  _et al._
(2018), we find an agreement between the observed and model melt ponds. To
match the melt pond statistics, we had to choose $r_{0}\approx 0.6\text{ m}$
in the model for both 1998 and 2005 ponds, in close agreement with $r_{0}$ we
find here from the height correlation functions of the pre-melt height
topographies (see Table 1). In all three measurements of the LiDAR topography,
we found $r_{0}$ within at most several cm from $r_{0}=0.6\text{ m}$. We take
this as evidence that melt ponds are controlled by the pre-melt snow
topography and that the “snow dune” model may be meaningfully extended to at
least the kilometer scale.

We will make several final notes to end this section. First, we note that
since $r_{0}\approx 0.6\text{ m}$ in observations for four different years and
multiple locations and ice types, the horizontal scale of the snow features
seems to be a robust property of the snow cover. For this reason, future
models that require a representation of snow may be able to keep the
correlation length $l_{0}\approx 5.6\text{ m}$ fixed, and only keep track of
mean snow depth and variance in order to represent snow in a principled way.
Second, we note that distributions other than a gamma distribution that have
been used to model the snow-depth, such as normal (Liston  _et al._ , 2018) or
lognormal (Landy  _et al._ , 2014), can also often fit the measurements well.
For example, a lognormal distribution passes the KS test for all three LiDAR
measurements, while a normal distribution passes the KS test for 2010
measurements. However, these distributions cannot qualitatively capture the
height distribution of the synthetic “snow dune” topography for all values of
the density of mounds, $\rho$, while a gamma distribution can. For small
$\rho$ (large ratio of standard deviation to mean), a normal distribution
predicts a significant fraction of negative observations, while a lognormal
distribution predicts a fat tail at large $h$. Since the “snow dune” model so
closely reproduces both the LiDAR and melt pond measurements, we believe that
a distribution that is qualitatively consistent with the predictions of the
“snow dune” model is a more justified model for a snow-depth distribution that
would fit the measurements over a wider parameter range. Resub.162Finally, we
note that Filhol and Sturm (2015) and Kochanski  _et al._ (2019a) observed an
approximately linear relationship between snow dune height and width for dunes
that formed on snow-covered lake ice, tundra, and mountain landscapes. Our
model prescription is therefore consistent with these observations. The data
in Kochanski  _et al._ (2019a) was quite noisy, but Filhol and Sturm (2015)
reported a slope of the height-width relationship ranging from about 0.03 to
0.06. We can estimate the mound aspect ratio in our model as $h_{m,0}/r_{0}$
and compare it with Filhol and Sturm (2015). From Table 1, we find that the
aspect ratio of mounds in our model varies between 0.01 for 2010 measurements
and 0.03 for 2009N measurements. These estimates are lower but broadly of the
same order of magnitude as those of Filhol and Sturm (2015), which suggests
that mounds in our model are at least broadly consistent with dunes observed
in the field, although not in detail.

##  V Heat transport through the ice

In this section we will consider the application of our model to heat
conduction through the ice. Sturm  _et al._ (2002) conducted measurements on
deformed multi-year ice, and concluded that the conductivity of snow inferred
from large-scale ice growth is roughly 2.4 times higher than the conductivity
of snow measured on-site. They ascribed a significant portion of this
difference to spatial variability of snow depth and also suggested that a
significant fraction of heat is transported horizontally. Motivated by their
study, here we investigate the extent to which these conclusions also apply to
undeformed ice. We consider a full 3-dimensional model of heat conduction to
determine how much heat is extracted through undeformed ice with variable snow
cover.

To model heat conduction, we assume that ice is a block of uniform thickness,
$H$, that ice and snow have fixed conductivities, $k_{i}$ and $k_{s}$, that
the snow cover is well-described by our “snow dune” model, that the
temperature at the ice-ocean interface is fixed at the freezing point of salt
water, $T_{f}$, that the temperature of the snow surface is fixed at the
temperature of the atmosphere, $T_{a}$, and that the temperature field within
the ice and snow is in a steady state. Additionally, we assume that
Resub.398heat is transported purely vertically in the snow (but not
necessarily in the ice), which is reasonable, since snow cover is typically
thin compared to its horizontal correlation length. Finally, we impose
periodic boundary conditions in the horizontal directions.
Resub.393Resub.M2.2.1Some of the assumptions above may not be realistic. For
example, ice is almost certainly not a uniform block, snow conductivity,
$k_{s}$, typically has large variations in both the horizontal and the
vertical, and the system is likely not in a steady state since the conditions
in the atmosphere change faster than snow and ice can adjust. For these
reasons, our investigation here should be considered to be a first order
estimate of the effects of snow geometry.

In SI section S3, we develop relations for the mean conductive heat flux
through the ice, $F_{c}$, under the assumptions above. In particular, starting
from Laplace’s equation for heat conduction and using the fact that our “snow
dune” model is fully characterized by the mean, variance and the correlation
length of the snow surface, we show that $F_{c}$ must be of the form

$\displaystyle F_{c}=F_{0}\Phi\big{(}\eta,\Sigma,\Lambda\big{)}\text{ , }$
(12) $\displaystyle F_{0}\equiv k_{i}\frac{T_{f}-T_{a}}{H}\text{ ,
}\eta\equiv\frac{k_{i}}{k_{s}}\frac{\langle h\rangle}{H}\text{ ,
}\Sigma\equiv\frac{\sigma(h)}{\langle h\rangle}\text{ ,
}\Lambda\equiv\frac{l_{0}}{H}\text{ , }$ (13)

where $\Phi$ is a non-dimensional flux that is determined by a non-dimensional
snow depth, $\eta$, a non-dimensional snow roughness, $\Sigma$, and a non-
dimensional correlation length, $\Lambda$. The dimensional quantity $F_{0}$ is
the heat flux conducted through the ice with no snow on top, so the function
$\Phi\leq 1$ represents the fraction of this flux that is conducted when snow
is present. After showing this, in SI section S3, we show that the heat
conduction problem can be explicitly solved if $\Lambda\rightarrow\infty$
(when heat transport is purely vertical), and when $\Lambda\rightarrow 0$
(when the horizontal heat transport dominates). Knowing these two limits, we
then numerically solve the heat conduction problem for various combinations of
the parameters $\eta$, $\Sigma$, and $\Lambda$ and show that the function
$\Phi$ is approximately equal to

$\displaystyle\Phi\big{(}\eta,\Sigma,\Lambda\big{)}\approx\Phi_{v}+\frac{\Phi_{h}-\Phi_{v}}{(1+c\Lambda)^{2}}\text{
, }$ (14)
$\displaystyle\Phi_{v}\equiv\int_{0}^{\infty}\frac{\tilde{f}_{\Gamma}(z,\Sigma)}{1+z\eta}\text{d}z\text{
, }$ (15) $\displaystyle\Phi_{h}\equiv\frac{1}{1+\eta(1-\Sigma^{2})}\text{ ,
}$ (16)

where $c\approx 0.83$ is a numerical constant, $\Phi_{v}$ is the non-
dimensional flux if vertical heat transport dominates, and $\Phi_{h}$ is the
non-dimensional flux if horizontal heat transport dominates. $\Phi_{h}$ is an
upper bound on heat flux given parameters $\eta$ and $\Sigma$. The gamma
distribution $\tilde{f}_{\Gamma}$ in Eq. 15 is normalized to have a mean equal
to 1, and is consequently only a function of $\Sigma$ with the parameters
$k=\Sigma^{-2}$ and $h_{0}=\Sigma^{2}$. Equation 16 is only valid for
$\Sigma<1$. If $\Sigma>1$, we have that $\Phi_{h}=1$ so that, if
$\Lambda\rightarrow 0$, the heat is conducted as if there were no snow on top
of the ice, $F=F_{0}$.

To quantify how much ice growth is due to snow-depth variability, we also
consider the non-dimensional flux assuming a uniform snow cover,
$\Phi_{u}\equiv\frac{1}{1+\eta}$. The fraction of ice growth that can be
attributed to snow-depth variability can then be estimated as
$\frac{\Phi-\Phi_{u}}{\Phi}$. We show the estimated values of all of the non-
dimensional parameters discussed above for the three LiDAR measurements we
considered in the previous section in Table 2.

From Eq. 14, we can see that the contribution of horizontal heat transport to
total heat transport is proportional to $(1+c\Lambda)^{-2}$. As we discussed
in the previous section, the correlation length for all datasets we considered
is around $l_{0}\approx 5.6\text{ m}$ and is likely on the same order
throughout the Arctic. This means that for ice of realistic thickness of
around $1\text{ m}$, $\Lambda\approx 5$. This implies that horizontal heat
transport contributes less than 5% of the maximum possible contribution,
$\Phi_{h}-\Phi_{v}$. This is very small - as we can see from Table 2, the
total non-dimensional flux, $\Phi$, is equal to $\Phi_{v}$ to two significant
digits for all measurements. Even for exceedingly thick ice of 2 m thickness,
and snow with significant variability, as in the 2009N measurements,
horizontal heat transport still contributes less than 1% of total ice growth.
Therefore, horizontal heat transport can most likely be ignored on flat
undeformed ice. This is in contrast with Sturm  _et al._ (2002) who found a
significant fraction of heat is transported horizontally in regions of
deformed ice. For this reason, it may be necessary to include non-uniform ice
thickness to appropriately capture heat transport on deformed ice.

Even though horizontal heat transport may be neglected on undeformed ice, snow
variability may still noticeably contribute to conductive heat flux. In Table
2, we show the fraction, $(\Phi-\Phi_{u})/\Phi$, of heat transport that is due
to snow-depth variability. We see that snow-depth variability contributes
between 4% for 2010 measurements that have the least snow-depth variability
and 11% for 2009N measurements that have the most snow-depth variability, with
the effect generally increasing with $\eta$ and $\Sigma$. Although this is not
large, it is comparable to the effect of melt ponds during summer - for
example, if 15% of ice of albedo 0.6 is covered by ponds of albedo 0.25, ice
will melt on average 10% faster than if there were no ponds. Therefore, the
effect of snow variability on first-year ice during winter may be enough to
partially compensate the effect of melt ponds during summer. Again, this
effect is smaller than on deformed ice, where Sturm  _et al._ (2002) found
that about a quarter of ice growth is due to variable snow and ice geometry.

Resub.M1.2We note that Eqs. 15 and 16, which quantify purely vertical and
purely horizontal heat flux, follow from general properties of Laplace’s
equation in the limits $\Lambda\rightarrow\infty$ and $\Lambda\rightarrow 0$,
and from the fact that the snow depth distribution is well-described by a
gamma distribution. These conditions do not change if we modify the “snow
dune” model details such as mound shape, anisotropy, or size distribution (SI
section S2.3), so our formula for total heat flux should be applicable
regardless of such specifics. The only factor that could potentially change
with these alterations is the detailed shape of the function that describes
the transition between vertical and horizontal heat transport regimes.

Resub.M2.2.2Resub.393BFinally, we note that if a reasonable statistical model
of spatial variability of snow conductivity can be developed, this variability
could be included in the above framework. Namely, if snow has a non-constant
conductivity, $k_{s}(\mathbf{x})$, heat transport would be affected by an
effective snow topography that accounts for these variations in conductivity,
$h(\mathbf{x})\frac{\langle k_{s}\rangle}{k_{s}(\mathbf{x})}$. Thus, a non-
constant snow conductivity would simply introduce an effective topography and,
with a good statistical model of $k_{s}(\mathbf{x})$, we could follow the same
steps to estimate the non-dimensional heat fluxes, $\Phi$, $\Phi_{v}$, and
$\Phi_{h}$. We note that if the horizontal length-scale of snow conductivity
variations is on the order of, or less than, ice thickness, it could induce
significant horizontal heat transport.

Table 2: Non-dimensional parameters of the snow surface for the LiDAR measurements of Polashenski  _et al._ (2012) assuming an ice thickness of $H=1\text{ m}$, conductivity of fresh ice $k_{i}=2.034\text{ Wm}^{-1}\text{K}^{-1}$ (Untersteiner, 1964), snow conductivity estimated by Sturm  _et al._ (2002), $k_{s}=0.14\text{ Wm}^{-1}\text{K}^{-1}$, and duration of stage I of pond evolution $T=5\text{ days}$.  | $\eta$ | $\Sigma$ | $\Lambda$ | $\Phi_{u}$ | $\Phi_{v}$ | $\Phi_{h}$ | $\Phi$ | $(\Phi-\Phi_{u})/\Phi$ | $\omega$
---|---|---|---|---|---|---|---|---|---
2009N | 2.2 | 0.5 | 5.5 | 0.31 | 0.35 | 0.38 | 0.35 | 0.11 | 1.7
2009S | 1.9 | 0.4 | 5.24 | 0.34 | 0.37 | 0.39 | 0.37 | 0.07 | 1.7
2010 | 1.9 | 0.3 | 5.75 | 0.34 | 0.36 | 0.37 | 0.36 | 0.04 | 1.5

##  VI Melt pond evolution during early melt season

During the early melt season, ice is largely impermeable and ponds can quickly
flood vast areas of the ice surface. This is known as stage I of pond
evolution and is typically followed by stages II and III which correspond to
pond drainage and subsequent slow pond growth on highly permeable ice
(Polashenski  _et al._ , 2012; Landy  _et al._ , 2014). In this section, we
will discuss equations for pond evolution during stage I when ice can be
considered to be impermeable, which we derive in SI section S4. These
equations follow from the fact that the pre-melt snow distribution on flat ice
is well-described by our “snow dune” topography.

###  VI.1 Analytic model for pond evolution during stage I

Since water can flow relatively freely through permeable snow, during stage I
of pond evolution there likely exists a common water table for the entire ice
floe, and ponds are the regions where surface topography lies below it.
Therefore, to model stage I of pond evolution, we will assume that snow is
fully permeable and that snow below the water level is completely saturated
with water, that underlying ice is impermeable and initially flat, that the
pre-melt snow topography is well-described by our “snow dune” model, and that
no meltwater is lost from the domain (we will reconsider this assumption in
the next subsection). The ice and snow topography can change by melting,
thereby creating meltwater. We will thus assume that snow, bare ice, ponded
ice, and A24ponded snow (water-covered snow) melt at different rates but that
these rates are constant in space and time and independent of factors such as
snow or pond depth. B20Finally, we will assume that both snow and ice melt
only at the surface and that ice cannot melt until snow has melted away.

Figure 5: a-d: Pond coverage evolution during stage I found using Eqs. 17 and
18 for a variety of model parameters. In each panel, we change one parameter,
while we keep the others fixed at default values $\dot{h}_{s}=4\text{ cm
day}^{-1}$, $r_{s}=0.4$, $\sigma(h)=0.05\text{ m}$, and $k=4$. a) Pond
coverage evolution for different snow melt rates, $\dot{h}_{s}$. b) Pond
coverage evolution for different snow to ice density ratios, $r_{s}$. c) Pond
coverage evolution on topographies with different snow-depth standard
deviations, $\sigma(h)$. d) Pond coverage evolution on topographies with
different shape parameters, $k$. e) Comparing pond evolution during stage I
calculated using Eqs. 17 to 19 and physical parameters described in the text
(dashed blue and yellow lines) with observations of Polashenski  _et al._
(2012) (red line). The dashed blue line represents pond evolution estimated
using Eqs. 17 and 18 that assumes no drainage during stage I. Dashed yellow
line represents pond evolution estimated using Eq. 19, which assumes a
constant drainage rate beyond the percolation threshold of $p_{c}=0.35$.
Topography used to calculate stage I evolution has the same mean and variance
as the actual measured pre-melt snow topography obtained with LiDAR
measurements.

In SI section S4, we develop equations for the time-evolution of pond
coverage, $p$, and water level, $w$, during stage I under the assumptions
above. We do this by tracking the volume of water produced by melting. We use
the fact that, if no water is lost from the domain, the water level can only
increase, so that ponds only encounter bare snow topography that melts at a
uniform rate and is, therefore, unaltered throughout the melt season. This
means that we can use our pre-melt topography to accurately describe the
ponded fraction of the surface. It also means that regions of bare ice (ice
with no snow and no water on top) do not exist. We greatly simplify the
problem by neglecting terms that are proportional to the density difference
between water and ice, which are, in addition to being small, only relevant at
very high pond coverage. Using this strategy, we show that pond coverage
evolution is approximately the solution to the system of ordinary differential
equations

$\displaystyle\dot{w}\approx-\frac{r_{i}r_{s}(1-p)}{1-r_{s}(1-p)}\dot{h}_{s}\text{
, }$ (17)
$\displaystyle\dot{p}=f_{\Gamma}(w-\dot{h}_{s}t)\Big{(}\dot{w}-\dot{h}_{s}\Big{)}\text{
, }$ (18)

where $\dot{x}\equiv\frac{\text{d}x}{\text{d}t}$ is a shorthand notation for
the rate of change of quantity $x$, $w$ is the water level, $p$ is the pond
coverage faction, $r_{i}\equiv\frac{\rho_{i}}{\rho_{w}}$ is the ratio of ice
density to meltwater density, $r_{s}\equiv\frac{\rho_{s}}{\rho_{i}}$ is the
ratio of snow density to pure ice density, $\dot{h}_{s}$ is the melt rate of
bare snow, assumed to be constant in space and time, and
$f_{\Gamma}(w-\dot{h}_{s}t)$ is the gamma distribution given by Eq. 9 and
evaluated at $w-\dot{h}_{s}t$. The parameters of the gamma distribution are
determined from the mean and variance of the pre-melt snow topography
according to Eq. 10. The initial conditions for Eqs. 17 and 18 are no ponds at
the initial time, $p(t=0)=0$, and a water level of zero, $w(t=0)=0$. Note that
$\dot{h}_{s}<0$, so $\dot{w}>0$ and $\dot{p}>0$. Even though Eqs. 17 and 18
are approximate, they are a very accurate approximation to the full 2d model
defined by the assumptions above at low pond coverage and only deviate
slightly from the full model at high pond coverage (see SI section S4). In SI
section S4, we also derive a second-order approximation that becomes nearly
indistinguishable from the full 2d model for any pond coverage.

We can see from Eqs. 17 and 18 that, since ice and water density are fixed,
stage I pond coverage evolution depends mainly on the density of snow through
$r_{s}$, the melt rate of bare snow, $\dot{h}_{s}$, and the mean and variance
of the initial snow-depth distribution through $f_{\Gamma}$. Note that melt
rates of bare ice, ponded ice, and ponded snow do not enter this first-order
approximation for pond evolution. In Figs. 5a-d, we change each of the
parameters that enter Eqs. 17 and 18 to see how they affect stage I pond
evolution. We note that in Figs. 5a-d, we parameterized $f_{\Gamma}$ with its
standard deviation, $\sigma(h)$, and its shape parameter, $k$, since
$\sigma(h)$ strongly affects pond coverage evolution, while $k$ affects it
only weakly. Within a reasonable range of parameters, pond evolution during
stage I is most sensitive to the rate of snow melt (Fig. 5a) and the standard
deviation of the initial snow-depth distribution (Fig. 5c). In fact, as we
explain in SI section S4, increasing the snow melt rate by some factor leads
to the same pond coverage evolution as decreasing the volume of snow by the
same factor. Snow density and the shape of the initial snow-depth
distribution, parameterized by $k$ and assuming that $\sigma(h)$ is fixed, do
not matter as much (Figs. 5b and d).

###  VI.2 Meltwater drainage during stage I

B21One important factor we have neglected is the potential for limited
drainage during stage I. Even though ice is typically impermeable during this
stage, outflow pathways such as cracks in the ice, seal breathing holes, or
the floe edge can exist (Polashenski  _et al._ , 2012; Fetterer and
Untersteiner, 1998; Eicken  _et al._ , 2002; Holt and Digby, 1985). Popović
_et al._ (2020) explain that if pond coverage is low, and the underlying
surface is impermeable, such isolated flaws in the ice should not
significantly affect the pond coverage. However, if pond coverage is above a
special value called the percolation threshold, $p_{c}$, there can be
significant pond drainage through these flaws. In fact, if the drainage rate
is great enough, the pond coverage would be limited to below $p_{c}$. For
example, unless ice deformation prevents water from flowing into the ocean,
the floe edge will limit pond growth to below $p_{c}$. Popović  _et al._
(2020) and Popović  _et al._ (2018) estimate $p_{c}$ to be between 0.3 and 0.4
for Arctic sea ice. Below, we show that including drainage is necessary to
accurately reproduce the evolution of pond coverage throughout stage I.

B21.2Including drainage in our model leads to some inconsistencies with the
assumptions under which we derived Eqs. 17 and 18. First, drainage opens the
possibility for the water level to decrease, $\dot{w}<0$, potentially exposing
bare ice and topography that was altered by different melt rates of different
regions of ice. These effects were not taken into account when deriving Eqs.
17 and 18, and, therefore, we can keep these equations only if the drainage is
not too great. In fact, Eqs. 17 and 18 remain valid when
$\dot{w}-\dot{h}_{s}>0$ (the altered topography remains submerged) and while
$w>0$ (there are no regions of bare ice). Second, the percolation threshold
only affects the pond coverage in the way we discussed above if the underlying
surface is impermeable. In our case, the surface is a combination of
impermeable ice and permeable snow. Since the pond coverage will likely exceed
the percolation threshold before all the snow melts, there likely exists a
complicated transition period between a surface that is a mix of permeable
snow and impermeable ice and a fully impermeable ice surface. Moreover, when
drainage is great enough so that the percolation threshold limits pond growth,
the water level in different ponds needs to be at least slightly different, in
contrast with the assumption of a common water table we made in the previous
subsection. Here, we simply ignore these complications and assume that Eqs. 17
and 18 are approximately valid even when pond drainage occurs and that
drainage is only activated if pond coverage exceeds the percolation threshold.
We leave a detailed analysis of pond drainage during stage I for another
study.

If we assume $\dot{w}-\dot{h}_{s}>0$ and $w>0$, we can straightforwardly
include drainage in pond evolution (Eqs. 17 and 18) by simply altering the
water balance equation. We derive this in SI section S4. For example, if ponds
drain at a rate of $Q$ cm per day, Eq. 18 remains the same and 17 is modified
to

$\dot{w}=-\frac{r_{i}r_{s}(1-p)\dot{h}_{s}+Q}{1-r_{s}(1-p)}\text{ . }$ (19)

To account for the fact that we expect little drainage to occur below the
percolation threshold, we can make $Q$ a function of $p$. The simplest
representation would be to include A26$Q(p)=Q_{0}\mathcal{H}(p-p_{c})$, where
$\mathcal{H}$ is the Heaviside function, equal to 0 or 1 depending on whether
its argument is less or greater than 0, and $Q_{0}$ is a constant that sets
the maximum drainage rate. Using this function yields the yellow line in Fig.
5e. If $\dot{w}-\dot{h}_{s}<0$, Eq. 18 predicts that the pond coverage will
decrease, $\dot{p}<0$. So, if $\dot{w}<\dot{h}_{s}$, we expect the pond
coverage to remain fixed at the percolation threshold, $p=p_{c}$.

In Fig. 5e we compare stage I pond evolution predicted by our model with the
measurements of stage I pond coverage evolution by Polashenski  _et al._
(2012). Since accurate measurements of pre-melt topography are available for
that time and location (Fig. 3), we can choose the parameters of the gamma
distribution highly accurately. The density of snow and the snow melt rate are
uncertain. We choose the density of snow to be $\rho_{s}=350\text{ kg
m}^{-3}$, and the albedo of melting snow to be $\alpha_{s}=0.74$, consistent
with observations of Polashenski  _et al._ (2012). We take the solar flux to
be $F_{\text{sol}}=254\text{ W m}^{-2}$, which Polashenski  _et al._ (2012)
report as the average flux of solar energy for the duration of the experiment.
We treat the sum of longwave, sensible, and latent heat fluxes, $F_{r}$, as a
tuning parameter. The value we choose, $F_{r}=-15\text{ W m}^{-2}$, is typical
of the region and season, and is therefore consistent with the measurements.
The blue line in Fig. 5e shows the stage I pond evolution predicted by Eqs. 17
and 18, while the red line represents the measurements of Polashenski  _et
al._ (2012). The two curves agree up to the point when pond coverage reaches
$p\approx 0.35$, a plausible value for the percolation threshold. Thus, the
measurements are consistent with drainage being activated beyond the
percolation threshold. A26.2To match pond evolution beyond $p=0.35$ we used
Eq. 19, and tuned the drainage rate, $Q_{0}$, to $2.6\text{ cm day}^{-1}$
(yellow line), B22a rate consistent with observations of Polashenski  _et al._
(2012). We can see that Eq. 19 is consistent with the observations, although
this matching cannot be fully confirmed due to uncertainty in the parameters.

###  VI.3 A condition for pond development

Figure 6: The horizontal axis shows the non-dimensional snow roughness,
$\sigma(h)/\langle h\rangle$, while the vertical axis shows the non-
dimensional water level, $\omega$, described in the text. The white-colored
region represents ice that remains pond-free throughout the summer, while the
dark-blue region represents ice that develops pond coverage that exceeds the
percolation threshold by the end of stage I. The light-blue region represents
ice that has some ponds by the end of stage I, but the pond coverage does not
exceed the percolation threshold.

Pond coverage reaches its peak by the end of stage I. Afterwards, ponds drain,
and, by the end of the drainage stage, the remaining ponds correspond to
regions of ice that are below sea level. Observations show that ponds that
remain after drainage are a subset of ponds that exist during peak pond
coverage (Polashenski  _et al._ , 2012; Landy  _et al._ , 2014). Therefore, if
ponds do not develop during stage I, they will likely never develop. Since the
pre-melt snow-depth distribution controls pond evolution during stage I, it,
therefore, also has a disproportionate effect on later pond evolution. If the
snow-depth is highly variable, ponds will develop sooner than if the snow
cover is uniform. If the snow cover is relatively uniform and stage I does not
last long enough, the ice may remain pond-free throughout the summer. Here, we
will show how we can use the model we developed above to derive a simple
criterion for whether ice will develop ponds or not.

We start with Eqs. 17 and 18. If we assume that pond coverage stays low
throughout stage I, we can set $p\approx 0$ in Eq. 17, and integrate it
directly to get the water level at time $t$ during stage I,
$w(t)\approx-\frac{r_{i}r_{s}}{1-r_{s}}\dot{h}_{s}t$. We can then non-
dimensionalize Eq. 18 by introducing non-dimensional parameters $\Sigma$ and
$\omega$ to get

$\displaystyle\frac{\text{d}p}{\text{d}\omega}=\tilde{f}_{\Gamma}(\omega,\Sigma)\text{
, }$ (20) $\displaystyle\Sigma\equiv\frac{\sigma(h)}{\langle
h\rangle}\quad\text{ , }$
$\displaystyle\quad\omega\equiv-\frac{1-r_{s}(1-r_{i})}{1-r_{s}}\frac{\dot{h}_{s}t}{\langle
h\rangle}\text{ , }$ (21)

Here, $\Sigma$ is the non-dimensional roughness as in section V and $\omega$
can be interpreted as a non-dimensional water level. As in Eq. 15 of section
V, $\tilde{f}_{\Gamma}(\omega,\Sigma)$, is a gamma distribution with a mean
equal to 1, so that it only depends on $\Sigma$, with parameters
$k=\Sigma^{-2}$ and $h_{0}=\Sigma^{2}$, and evaluated at $\omega$. Thus, from
Eq. 20, we get a simple relation

$p=F_{\Gamma}(\omega,\Sigma)\text{ , }$ (22)

where $F_{\Gamma}(\omega,\Sigma)$ is a cumulative gamma distribution with a
mean equal to 1,
$F_{\Gamma}(\omega,\Sigma)\equiv\int_{0}^{\omega}\tilde{f}_{\Gamma}(z,\Sigma)\text{d}z$.
Therefore, if pond coverage is low, it depends only on $\omega$ and $\Sigma$
and we can express it simply in terms of the cumulative gamma distribution.

Let us now assume that stage I lasts for some time, $T$, on the order of 5
days. If the ponds do not develop within that time, the ice will remain pond-
free throughout the summer. We can set some threshold, $p^{*}$, say
$p^{*}=0.01$, and consider ice to be pond-free if the pond coverage is below
$p^{*}$. We can then find the boundary in $\omega$-$\Sigma$ space,
$\omega^{*}(\Sigma)$, such that if $\omega<\omega^{*}$, ice remains pond-free
at the end of stage I. According to Eq. 22, this boundary is given by

$\omega^{*}(\Sigma)=F^{-1}_{\Gamma}(p^{*},\Sigma)\text{ , }$ (23)

where $F^{-1}_{\Gamma}$ is an inverse of the cumulative gamma distribution
(gamma percentile function) with a mean equal to 1, $k=\Sigma^{-2}$, and
$h_{0}=\Sigma^{2}$, evaluated at $p^{*}$. Equation 23 gives a universal
criterion that ice remains pond-free. If $\omega$, calculated given the snow
melt rate, density, depth, roughness, and the duration of stage I, exceeds
$\omega^{*}$, ice will develop at least some ponds. The boundary,
$\omega^{*}(\Sigma)$, is largely insensitive to the choice of the threshold
coverage, $p^{*}$, so long as $p^{*}$ is much smaller than 1.

In the case when ice develops some ponds, we can also ask how developed those
ponds will be. In particular, we can assume that ponds are fully developed if
the pond coverage exceeds the percolation threshold by the end of stage I.
Beyond the percolation threshold, pond drainage likely starts to play an
important role, and the post-stage I evolution of fully-developed ponds likely
proceeds in a fairly typical manner. Therefore, there likely exist three
typical trajectories for pond evolution during summer - pond-free, fully-
developed, and intermediate. To see whether the ponds will develop fully, we
can again use Eq. 23 but use $p^{*}=p_{c}\approx 0.35$. Using Eq. 23 with
$p^{*}=0.35$ is approximate since the pond coverage is no longer small, but it
is still reasonably accurate.

In Fig. 6, we show the conditions for these three regimes of pond development.
We can see that when $\omega>1$, ponds always fully develop since, in this
case, the water level by the end of stage I exceeds the height of snow even
when the snow cover is uniform. Otherwise, if $\omega<1$, ponds develop more
easily when the snow-depth has more variability. We can estimate $\omega$ and
$\Sigma$ for the measurements of Polashenski  _et al._ (2012). We use
$T=5\text{ days}$, snow depth mean and variance estimated from the LiDAR
measurements, and the same snow melt rate density as in the previous
subsection. We show $\Sigma$ and $\omega$ estimated in this way in Table 2. We
can see that $\omega>1$ in all cases so that we expect that ponds develop
fully. In the study of Polashenski  _et al._ (2012) this indeed happened. In
order for ice to have remained pond-free, stage I would have had to have
lasted less than 1.5 days in 2010, which had the most uniform snow, and less
than 0.5 days for 2009N, which had the most variable snow. We note that
persistently pond-free ice has been observed in the Arctic (Perovich  _et al._
, 2002), while Antarctic sea ice only rarely develops ponds. Resub.655Although
we do not have snow measurements from these pond-free ice floes, general
conditions on Antarctic ice are broadly consistent with our framework here.
Snow on Antarctic ice is typically deeper than in the Arctic, lowering both
$\Sigma$ and $\omega$, and it is both thinner and more saline than Arctic ice,
potentially lowering the duration of stage I, $T$, thereby further lowering
$\omega$. It would be interesting to see where these pond-free ice floes fall
within our pond-development diagram.

##  VII Discussion

Resub.656Resub.M.3Throughout this paper, we showed that our model can
reproduce key statistical features on three nearby locations where detailed
LiDAR data were collected as well as two other locations where melt pond
images were taken. Resub.665Although this is not enough to claim generality,
our results are encouraging, suggesting that our model may be applicable in
large-scale studies of snow on undeformed Arctic sea ice and could serve as a
starting point for studying snow on other terrains such as deformed ice,
Antarctic ice, or Alpine snow fields. Resub.667Moreover, as the fraction of
first-year ice increases in the Arctic due to global warming, the conditions
similar to the ones we tested here will likely become more common. Here, we
will briefly discuss how our results can be applied to large-scale studies,
assuming their generality holds.

  1. 1.

We believe that large-scale models should use a gamma distribution for snow
depth on undeformed ice. Even though other distributions, such as normal or
lognormal, may still work well enough, adopting a gamma distribution is a more
principled approach that likely works over a wider range of parameters.
Therefore mean and variance of snow depth are enough to characterize the whole
distribution according to Eqs. 10. Resub.M.4Since, according to the LiDAR as
well as melt pond data, the correlation length appears to be stable under a
variety of environmental conditions, it may be possible to use our “snow dune”
model with a horizontal scale parameter $r_{0}=0.6\text{ m}$ ($l_{0}=5.6\text{
m}$) to study horizontal snow distribution.

  2. 2.

When modeling heat conduction through undeformed ice, it is likely safe to
neglect horizontal heat transport. On deformed ice, the situation may be more
complicated and further study is required.

  3. 3.

To parameterize ice albedo during summer, Eqs. 17 and 18 may be used for pond
evolution during stage I. Without knowing the outflow rate, $Q$, a reasonable
approach would be to cap pond coverage to below the percolation threshold of
0.3 to 0.4. Equations 17 and 18 can be combined with similar approaches, such
as the model of Popović  _et al._ (2020) for stages II and III, or the model
of Popović  _et al._ (2018) for stage III, to yield an analytic and accurate
model of pond evolution throughout the entire melt season.

  4. 4.

Resub.676Eq. 23 provides a yet-untested prediction of whether ice will remain
pond-free throughout the summer - we were not able to test this criterion
against real data appropriately, as in all three cases we considered, the
ponds have fully developed. So, it would be interesting to see whether this
criterion can explain the observations of pond-free ice in the Arctic and the
Antarctic. The non-dimensional water level, $\omega$, in Eq. 23 depends on the
duration of stage I. So, to be able to use this equation in large scale-
studies, it will be necessary to relate the duration of stage I to parameters
that are available on the large scale, such as energy fluxes, thickness, or
the salinity of ice.

##  VIII Conclusions

Snow cover greatly impacts sea ice evolution. It insulates the ice during the
winter, reflects sunlight during the summer, and controls melt pond evolution.
Resub.685In this paper, we presented an idealized geometric model of “snow
dune” topography that was capable of capturing key statistical properties of
the snow cover on undeformed first-year sea ice on three sites for which
detailed LiDAR measurements were available as well as for two other missions
that collected airborne melt pond images. The main conclusion of our paper is
that, at least for the measurements we considered, our “snow dune” model is a
very accurate representation of the 3-dimensional snow distribution on flat,
undeformed Arctic sea ice.

By comparing the moments of the height and snow-depth distributions, we showed
that the height distribution of our model is statistically indistinguishable
from the snow-depth distribution measured in detailed LiDAR scans on flat
undeformed ice. We also compared the correlation functions for the model and
the measurements to show that the horizontal statistics of snow are also well-
captured by our model. In addition to this, our model captures melt pond
geometry derived from helicopter images that span areas of 1 km2 with the same
model parameters as the LiDAR scans. This suggests that our model may capture
the horizontal properties of snow on such large scales. We note that in all
cases we tested, the correlation length of horizontal features (either snow or
melt pond), was around $l_{0}\approx 5.6\text{ m}$, suggesting that this
property Resub.701may be constrained throughout the Arctic. Since our model is
fully defined by only three parameters, it follows that it is enough to know
the mean snow depth, its variance, and the horizontal correlation length to
fully characterize the statistics of snow on undeformed ice. On ice that was
slightly deformed we showed that the measured snow-depth distribution contains
subtle differences from the model and that there exist long-range correlations
in the horizontal, inconsistent with our model.

After comparing our “snow dune” model to measurements, we considered its
application to heat conduction through the ice and to the development of melt
ponds. We solved a 3-dimensional heat conduction equation assuming the ice is
flat and snow cover that is well-described by our model. We developed simple
formulas that solve the conductive heat flux problem for any configuration of
snow consistent with our model. Using these formulas, we showed that,
Resub.713for any realistic ice thickness and snow dune size, horizontal heat
transport may be neglected on flat, undeformed ice (although this may not be
the case for deformed ice). We then also used our model to develop a simple
model for melt pond evolution during early stages of pond development when ice
can be considered impermeable. We find that pond coverage evolution can be
estimated using a system of two coupled ordinary differential equations that
closely approximate an equivalent 2d model, but are much simpler to solve.
They reveal connections with measurable environmental parameters. For example,
we show that doubling the amount of solar radiation has approximately the same
effect on pond growth rate as halving the snow depth. We also show that they
are consistent with observations. Using these equations, we develop a
criterion in terms of non-dimensional properties of the ice surface for
whether ice will develop ponds or remain pond-free during summer.

In this paper, we demonstrated that Resub.725many of the details of ice
surface conditions can be summarized by only a few parameters if one is
interested in bulk properties such as the total heat conducted or the mean
melt pond coverage. Resub.728Of the three parameters that define our model, it
may be the case that only two change from situation to situation, as the
correlation length may be stable, as discussed above. This means that future
GCMs may need to only keep track of the total amount of snow and the snow
depth variance in order to be able to represent the snow surface in a
principled way. Our results, therefore, uncover an important property of sea
ice snow cover that may improve the realism of sea ice models Resub.733while
adding little to their complexity.

###### Acknowledgements.

We thank BB Cael and Stefano Allesina for reading the paper and giving
comments. We thank Don Perovich for providing the melt pond image in Fig. 1d.
Predrag Popović was supported by a NASA Earth and Space Science Fellowship.
This work was partially supported by the National Science Foundation under NSF
award number 1623064 and by National Science Foundation award DMS-1517416
(MS). Code to generate the synthetic “snow dune” topography is archived as
Popović (2020) and is also available at https://github.com/PedjaPopovic/Snow-
dune-topography. LiDAR measurements collected by Polashenski  _et al._ (2012)
are available for download at http://chrispolashenski.com/data.php or
https://data.eol.ucar.edu/dataset/106.426. The authors declare no conflict of
interest.

## References

  * Popović  _et al._ (2018) P. Popović, B. Cael, M. Silber, and D. S. Abbot, Physical review letters 120, 148701 (2018).
  * Perovich and Richter-Menge (2009) D. K. Perovich and J. A. Richter-Menge, Annual Review of Marine Science 1, 417 (2009).
  * Stroeve  _et al._ (2007) J. Stroeve, M. M. Holland, W. Meier, T. Scambos, and M. Serreze, Geophysical research letters 34 (2007).
  * Holland and Curry (1999) M. M. Holland and J. A. Curry, Journal of climate 12, 3319 (1999).
  * Perovich (1996) D. K. Perovich,  _The optical properties of sea ice._ , Tech. Rep. (Cold Regions Research and Engineering Lab Hanover NH, 1996).
  * Yen (1981) Y.-C. Yen,  _Review of thermal properties of snow, ice and sea ice_ , Tech. Rep. (Cold Regions Research and Engineering Lab Hanover NH, 1981).
  * Sturm  _et al._ (2002) M. Sturm, D. K. Perovich, and J. Holmgren, Journal of Geophysical Research: Oceans 107, SHE (2002).
  * Polashenski  _et al._ (2017) C. Polashenski, K. M. Golden, D. K. Perovich, E. Skyllingstad, A. Arnsten, C. Stwertka, and N. Wright, Journal of Geophysical Research: Oceans 122, 413 (2017).
  * Mather (1962) K. Mather, Polar Record 11, 158 (1962).
  * Filhol and Sturm (2015) S. Filhol and M. Sturm, Journal of Geophysical Research: Earth Surface 120, 1645 (2015).
  * Petrich  _et al._ (2012) C. Petrich, H. Eicken, C. M. Polashenski, M. Sturm, J. P. Harbeck, D. K. Perovich, and D. C. Finnegan, Journal of Geophysical Research: Oceans 117 (2012).
  * Liston  _et al._ (2018) G. E. Liston, C. Polashenski, A. Rösel, P. Itkin, J. King, I. Merkouriadi, and J. Haapala, Journal of Geophysical Research: Oceans 123, 3786 (2018).
  * Kochanski  _et al._ (2019a) K. Kochanski, R. S. Anderson, and G. E. Tucker, The Cryosphere 13, 1267 (2019a).
  * Sharma  _et al._ (2019) V. Sharma, L. Braud, and M. Lehning, The Cryosphere 13, 3239 (2019).
  * Kochanski  _et al._ (2019b) K. Kochanski, G.-C. Defazio, E. Green, R. Barnes, C. Downie, A. Rubin, and B. Rountree, Journal of Open Source Software 4, 1699 (2019b).
  * Polashenski  _et al._ (2012) C. Polashenski, D. Perovich, and Z. Courville, Journal of Geophysical Research: Oceans 117 (2012).
  * Landy  _et al._ (2014) J. Landy, J. Ehn, M. Shields, and D. Barber, Journal of Geophysical Research: Oceans 119, 3054 (2014).
  * Untersteiner (1964) N. Untersteiner, Journal of Geophysical Research 69, 4755 (1964).
  * Fetterer and Untersteiner (1998) F. Fetterer and N. Untersteiner, Journal of Geophysical Research: Oceans 103, 24821 (1998).
  * Eicken  _et al._ (2002) H. Eicken, H. Krouse, D. Kadko, and D. Perovich, Journal of Geophysical Research: Oceans 107, SHE (2002).
  * Holt and Digby (1985) B. Holt and S. A. Digby, Journal of Geophysical Research: Oceans 90, 5045 (1985).
  * Popović  _et al._ (2020) P. Popović, M. C. Silber, and D. S. Abbot, Journal of Geophysical Research: Oceans 125, e2019JC016029 (2020).
  * Perovich  _et al._ (2002) D. Perovich, W. Tucker, and K. Ligett, Journal of Geophysical Research: Oceans 107 (2002).
  * Popović (2020) P. Popović, “PedjaPopovic/Snow-dune-topography: Code to create an artificial ’snow dune’ topograhy,”  (2020).