# Supporting Information for “Critical percolation threshold restricts late-
summer Arctic sea ice melt pond coverage”

††journal: Journal of Geophysical Research: Oceans

Introduction

In this document, we first discuss details needed to reproduce the model
results and derive equations stated in the main text and then provide some
additional results and discussions. This document is organized as follows. In
section S1, we describe the numerical implementation of the hole model; in
section S2, we discuss model assumptions; in section S3, we describe different
models used to generate synthetic topographies; in section S4, we describe the
relationship between the parameters in Eqs. 12 and 14 of the main text and
parameters of the full 2d model; in section S5, we describe the measurements
used to estimate model parameters; in section S6, we derive the universal
function, $g(\eta)$; in section S7, we derive the hole opening timescale,
$T_{\text{h}}$; in section S8, we show how the minimum pond coverage,
$p_{\text{min}}$, depends on the measurable properties of sea ice; finally, in
section S9, we discuss the hole opening distribution, $F$, and show that only
its tail affects the pond evolution.

##  S1 Details of the hole model

We implemented our model on a grid, typically $500\times 500$ grid points in
size. For each grid point, $\mathbf{x}$, we tracked two fields - ice surface
height, $h(\mathbf{x})$, and water level, $w(\mathbf{x})$. Points where
$w(\mathbf{x})>h(\mathbf{x})$ were defined as ponds. We assumed that each
connected pond has the same water level, but different ponds can have
different water levels. For convenience, we defined $w$ to be equal to ice
height for non-ponded grid points, $w(\mathbf{x})=h(\mathbf{x})$. We define
the sea level to be the origin of the vertical axis, so that $w(\mathbf{x})=0$
signifies that ponds are at sea level. We initialized the model by generating
a random surface, $h(\mathbf{x})$, according to a topographic model (see
section S3), and setting an initial water level. Typically, we completely
covered the surface with water by prescribing a water level equal to the
maximum ice height. In this way, the pond coverage fraction at the beginning
of a typical simulation was equal to 1. We also prescribed an initial ice
thickness, $H$, assumed to be uniform in our model, and shifted $w$ and $h$ by
an equal and constant amount to enforce hydrostatic balance. For a water level
field as defined above, with $w(\mathbf{x})=h(\mathbf{x})$ for non-ponded
regions, the hydrostatic balance constraint is simply

$\langle w(\mathbf{x})\rangle=\frac{\rho_{w}-\rho_{i}}{\rho_{w}}H\text{ , }$
(S1)

where $\langle...\rangle$ stands for the average over all grid points, and
$\rho_{w}$ and $\rho_{i}$ are water and ice densities. Before running the
model, we scaled the ice surface height field to have a standard deviation of
at most 2% of ice thickness, $H$. This ensured that suddenly removing all the
water from the ice surface would leave no depressions below sea level after
hydrostatic balance is enforced.

We considered each grid point to be a potential hole, and to each one, we
ascribed a “critical temperature,” $\theta_{c}(\mathbf{x})$, above which it
opens. These critical temperatures were independently drawn from a prescribed
probability distribution, $f_{\theta}$. As a default, we used a normal
$f_{\theta}$. The mean and variance of this distribution are arbitrary, so we
set the mean to 0 and variance to 1. We kept track of the “bulk ice
temperature,” $\theta$. All the grid points that had a critical temperature
below the current bulk ice temperature were considered open. We set the
initial ice temperature to be below the critical temperature of all grid
points. We note that the “critical temperature,” $\theta_{c}(\mathbf{x})$ and
the “bulk ice temperature,” $\theta$, are simply a way to diagnose where the
holes open in our model and need not be related to actual ice temperature.

A model step consisted of several sub-steps:

  1. 1.

First, we increase the bulk ice temperature, $\theta$, by a fixed amount,
$\text{d}\theta$, and open holes at locations where the new bulk temperature
exceeds the critical temperature. Multiple holes can open at each step. The
amount, $\text{d}\theta$, by which temperature increases within a given
interval of time, $\text{d}t$, depends on the pre-defined hole opening
timescale, $T_{\text{h}}$. This timescale is defined such that the bulk
temperature increases by one standard deviation of the prescribed critical
temperature distribution, $f_{\theta}$, during $T_{\text{h}}$. Since we set
the standard deviation of $f_{\theta}$ to 1, this means that we choose
$\text{d}\theta=\frac{\text{d}t}{T_{\text{h}}}$.

  2. 2.

Next, we drain the ponds through these newly opened holes. We simulate the
drainage iteratively. We first check which holes are active, i.e., are ponded
($w(\mathbf{x})>h(\mathbf{x})$) and are above sea level ($w(\mathbf{x})>0$).
Next, we find the ponds connected to these active holes. Then, we decrease the
water level for each of these ponds by a small amount. This changes the pond
coverage slightly, and we find the new ponds as regions with
$w(\mathbf{x})>h(\mathbf{x})$, and enforce $w(\mathbf{x})=h(\mathbf{x})$
otherwise. Again, we identify holes that are still active, find ponds
connected to these holes, and repeat the drainage process until all the holes
become inactive. A hole becomes inactive when either a hole’s pond water level
reaches the ice height at the location of the hole,
$w(\mathbf{x})=h(\mathbf{x})$, or when pond water level reaches sea level,
$w(\mathbf{x})=0$.

  3. 3.

Once drainage through all of the newly opened holes finishes, we
preferentially melt ponded ice. We do this by decreasing ice surface height,
$h$, in ponded regions while leaving it unchanged in bare ice regions. The
amount by which the ponded ice melts depends on the prescribed melt rate,
$\text{d}h_{\text{diff}}/\text{d}t$. We can neglect the melting of bare ice
since it is only the relative melt that changes the topography and affects
pond evolution. This is justified because we treat ice thinning independently
from surface melt, by prescribing a separate parameter, $\text{d}H/\text{d}t$,
that we can change irrespective of $\text{d}h_{\text{diff}}/\text{d}t$, as we
explain below.

  4. 4.

Finally, we enforce Eq. S1 to maintain hydrostatic balance. We do this by
shifting the entire ice surface by some amount. Water level is also shifted by
the same amount except for ponds at sea level that contain a hole and are
maintained at sea level. By returning ice to hydrostatic balance, we
effectively adjust the sea level relative to mean ice height.

After adjusting the hydrostatic balance a model step is complete. Then the
whole process repeats until all the ponds are at sea level. In the original
version of the model, we neglected ice thinning, but when studying stage III
of pond evolution, we also thinned the ice using a prescribed rate,
$\text{d}H/\text{d}t$, before enforcing hydrostatic balance. In the
simulations that include ice thinning, we treat the ice thinning rate as an
independent parameter from the surface melt rate, so that
$\text{d}h_{\text{diff}}/\text{d}t$ controls topography evolution while
$\text{d}H/\text{d}t$ controls ice thickness evolution, and we do not
prescribe pond and bare ice melt rates individually. We also tested versions
of the model where not every grid point contains a hole, and where more than
one hole can form at each grid point. These modifications did not affect our
conclusions in any noticeable way.

In summary, Table S1 shows variables and parameters that enter our model.

Table S1: A table of variables and parameters of the hole model. Variable |  Name
---|---
$\mathbf{x}$ |  Coordinate
$h(\mathbf{x})$ |  Ice height above sea level
$w(\mathbf{x})$ |  Water level
$\theta$ |  Bulk ice temperature
$\theta_{c}(\mathbf{x})$ |  Critical temperature of individual brine channels
$f_{\theta}$ |  Critical temperature distribution
$\text{d}h_{\text{diff}}/\text{d}t$ |  Melt rate of ponded ice relative to bare ice
$T_{\text{h}}$ |  Hole opening timescale
$H$ |  Ice thickness
$\text{d}H/\text{d}t$ |  Ice thinning rate

##  S2 Model assumptions

Our model was aimed at explaining the observation that ponds organize around
the percolation threshold and exploring the implications of that observation.
As such, the model contained only the necessary elements to explain the
organization around the percolation threshold, and we have made several
simplifying assumptions. All of these assumptions are likely well-justified
for undeformed first-year ice, and our model seems to be a good first-order
representation of pond evolution after the beginning of drainage. However, in
some cases the assumptions we made may not hold. We now discuss violations of
our assumptions and how some of them might be accounted for in a future
version of the model.

  1. 1.

First, by separating pond drainage from the melting and hole opening steps, we
have effectively assumed that drainage happens instantaneously. This
assumption is supported by observations that find that large ponds can be
drained in a matter of hours [11], whereas melting through the thickness of
the freeboard occurs on the timescale of several days. One reason why this
assumption may not hold is that it takes time for brine channels to enlarge to
macroscopic sizes. Although there are no measurements of the hole enlargement
process, a simple model by Polashenski et al. [11] suggests that this may
occur on timescales of roughly a day, comparable to the timescale of melt.
Another reason is that ponds are partially replenished by additional meltwater
flowing into them. The rate of drainage depends primarily on the hole radius,
the hydraulic head and the number of holes draining a particular pond, while
the rate of replenishing scales with the ice melt rate and the area of the
pond. So, large ponds replenish water more easily and it may be necessary to
wait until several holes open and grow to a certain size before such ponds can
actually start draining. For this reason, our assumption of quick drainage may
be invalid especially at the beginning of stage II when ponds cover a large
area. We tested our model including a non-negligible time to drain ponds,
$T_{\text{d}}$. In this case we did not separate the drainage and melting
steps, but rather melted the ice and adjusted the hydrostatic balance after
each increment of drainage, and we assumed that each pond requires
$T_{\text{d}}$ time to drain. Our simulations included drainage timescales of
up to 60% of the hole opening timescale, $T_{\text{h}}$, and up to 30% of the
memorization timescale, $T_{\text{m}}$. We found that simulations with long
drainage timescales deviated somewhat from the universal function, $g(\eta)$,
with the post-drainage pond coverage deviating up to $10\%$ from the
instantaneous drainage predictions. Therefore, although this effect is likely
not of primary importance, it may not be negligible during the initial
drainage period.

  2. 2.

Second, we assumed that ice is very flat at the beginning of stage II. The
topography of undeformed first-year ice that underlies the snow cover at the
beginning of stage I is likely very flat. However, it is questionable whether
first year ice at the beginning of stage II can also be considered to be flat
since different rates of melt of ponded ice, bare ice, and snow during stage I
may act to amplify topography variations. Non-negligible variations in the
topography will change the time it takes for different regions of ice to fall
below sea level. For this reason, it may be necessary to include the effects
of non-negligible initial topography variance in order to get accurate
estimates of pond evolution.

  3. 3.

Third, we assumed that none of the physical parameters change with time, which
may not be correct. In particular, the two timescales, $T_{\text{m}}$, and,
$T_{\text{h}}$, both depend on parameters such as the solar flux or ice
thickness that are not constant. Furthermore, the warming rate of the ice
interior that enters the hole opening timescale depends on the temperature of
the ice, and so changes as the ice warms. Short-time fluctuations, such as
variation due to daily changes in solar flux, likely do not contribute
significantly to pond evolution. However, longer time variation in these
parameters likely affect pond evolution quantitatively, although qualitative
conclusions likely remain the same. Time variability can be simply added to
our model. In this case, instead of using Eqs. LABEL:eqn:12 and LABEL:eqn:14
of the main text, pond coverage evolution would be found by solving a
differential equation

$\frac{\text{d}\Pi}{\text{d}t}=\frac{\text{d}g}{\text{d}\eta}\Big{(}\frac{\partial\eta}{\partial
t}+\sum_{i}\frac{\partial\eta}{\partial P_{i}}\frac{\partial P_{i}}{\partial
t}\Big{)}\text{ , }$ (S2)

where $P_{i}$ stand for all the parameters that change with time. The term
$\text{d}g/\text{d}\eta$ is the derivative of the universal function with
respect to its argument and is approximately given by Eq. LABEL:eqn:4 of the
main text.

  4. 4.

Next, we assumed that water flowing towards the holes does not change the ice
topography. Faster flowing water will exchange heat with the underlying ice at
a higher rate. This will then lead to the formation of channels near the
drainage holes that focus water and enhance melting. These channels have been
observed in the field near large holes [11, 8]. They have the potential to
impact the pond evolution by allowing a single hole to drain larger portions
of the ice. To assess the effect of channel formation, we added a crude
channel formation scheme to our model. To allow for channel formation, we
increased the melt rate along the medial line of ponds while they were
draining. To include some elements of realism, we made this increase in melt
rate inversely proportional to the distance from the hole, inversely
proportional to the local width of the pond in the transverse direction of the
medial line, and proportional to the volume of water flowing through the hole,
assuming that the water volume is equally partitioned between each hole that
drains a particular pond. Such a setup led to the formation of channels
between ponds with the most intense channeling occurring during initial
drainage since the highest water volume is drained at that time. These
assumptions are not highly realistic but we believe they are sufficient for
this initial test. We ran several such simulations changing the intensity of
melt along the medial line. We found that for a high enough channeling rate, a
single hole may be enough to drain the entire surface. However, for any
intermediate value, we found that the curve $g(\eta)$ remains approximately
the same, and the effect of channel formation can be summarized by simply
adjusting the values of $p_{c}$ and $c$. In particular, the highest channeling
rate we tested that did not immediately drain the whole surface reduced
$p_{c}$ on a “diffusion” topography from $p_{c}=0.5$ to $p_{c}\approx 0.27$
and $c$ from $c\approx 4.1$ to $c\approx 1.4$. This suggests that Eqs.
LABEL:eqn:12 and LABEL:eqn:14 of the main text may still be used in the
presence of channel formation, but the values of $p_{c}$ and $c$ may depend on
channel formation physics. The fact that the effect of channel formation can
be summarized by simply adjusting the values of $p_{c}$ and $c$ in our scheme
is likely due to the fact that most channeling occurs after drainage through
the first several holes, so the topography remains relatively unchanged as
more holes open. This likely also occurs in the field, but in a scenario where
the surface is modified throughout pond drainage, the values of $p_{c}$ and
$c$ may need to be made time-dependent as well. Exploring this hypothesis is
beyond the scope of this paper. We note that Popović et al. [14] estimated
that the percolation threshold for late-summer ponds photographed during the
1998 SHEBA mission was around 0.3, while it was around 0.4 for ponds
photographed during the 2005 HOTRAX mission. This is in contrast with the
“snow dune” topography which accurately describes the pre-melt surface
conditions on first-year ice [15] and predicts a percolation threshold of
roughly between 0.4 and 0.5, with values between 0.45 and 0.5 being more
likely based on the LiDAR measurements of the pre-melt snow topography [15].
Therefore, it is possible that channel formation by water flow decreases the
percolation threshold of the pre-melt snow topography by an amount on the
order of 0.1 to 0.2.

  5. 5.

We also assumed that at the beginning of stage II, ice is completely covered
by water. This assumption is unrealistic as ice is typically covered
significantly less at peak coverage. However, this assumption does not change
any of the conclusions we made about pond evolution. As long as the pond
coverage fraction at the beginning of the drainage stage is above the
percolation threshold, pond evolution will progress identically after drainage
below the percolation threshold following the opening of the first several
holes. The only effect of this assumption is that non-uniform melt rates
before the beginning of drainage may amplify the topography variations as
already discussed in item 2 above.

  6. 6.

Next, we assumed that drainage happens entirely through large holes and not
through the bulk of the ice. This is in contrast with all previous melt pond
models that modeled melt ponds as a balance between meltwater production and
drainage through the bulk of the ice [see, e.g., 9, 18]. Polashenski et al.
[11] observed that during stage II, water is almost entirely drained through
large holes. Furthermore, Eicken et al. [5] found that early in the melt
season, ice permeability is negligible. However, Polashenski et al. [12] found
that very large brine channels may not close fully, and there may be some
leftover bulk permeability of the ice. A significant amount of drainage
through the bulk of the ice during stage II could in principle qualitatively
alter the conclusions of our model. In particular, since the percolation
threshold does not control bulk drainage, $p_{c}$ would no longer be a
significant coverage fraction and the universal curve $g(\eta)$ would no
longer control pond evolution. However, since direct observations by Popović
et al. [14] show that ponds seem to be organized around the percolation
threshold, bulk ice drainage likely has only a small effect.

  7. 7.

We assumed that ice is a rigid plate that cannot elastically flex as the ponds
drain. This implies that hydrostatic balance is a global condition that
determines the height of the freeboard, rather than being determined by local
mass-balance. Equation LABEL:eqn:5 of the main text follows from this
assumption. This assumption is likely well-justified as ice flexure following
pond drainage has not been observed in field. The rigid-plate approximation is
also employed in all current melt pond models.

  8. 8.

When simulating stage III, we assumed that ponds grow entirely due to ice
thinning and we have neglected lateral melting. This assumption seems
justified at least on first-year ice, as observations show that freeboard
sinking due to ice thinning is the dominant mode of pond evolution during
stage III [11, 8], although some models suggest lateral melting may be
important [17]. Lateral melt would increase pond coverage during stage III,
and could be included in our model by adding a lateral melt contribution as
described in Popović and Abbot [13].

  9. 9.

We assumed that ice represents a rigid barrier for water flowing horizontally.
This is generally true, although water may be able to flow between
disconnected ponds if they are only separated by a narrow strip of ice. This
was demonstrated by Eicken et al. [4] who showed that tracers released in
disconnected ponds can actually get mixed. However, this effect is most likely
small. We note that while ice is mainly impermeable to horizontal flow of
water, snow is not, so our “snow dune” topography is actually meant to
represent ice topography left after snow has melted away.

##  S3 Synthetic topographies

To test our model, we used several types of synthetic topographies. Each of
these topographies had different statistical characteristics such as the
surface height distribution and the percolation threshold which allowed us to
test the generality of our results. Importantly, each of these topographies
was well-described by a single characteristic length-scale, which is necessary
in order for the universal function $g$ to describe the drainage well.

Our default was the “snow dune” topography (Fig. S1a). Ponds that form on this
surface are shown in Fig. S1d. This topography is a generalization of the void
model described in Popović et al. [14], and ponds that form on this topography
therefore reproduce the pond geometry well. In Popović et al. [15], we also
showed that this topography reproduces the measured snow topography on
undeformed ice highly accurately. We note that even though this topography was
meant to describe the properties of the pre-melt snow surface, here we will
consider it as a representation of impermeable ice. We generate this
topography as a sum of $N_{\text{m}}$ mounds placed randomly on an initially
flat surface. These mounds have a Gaussian shape,
$h(\mathbf{x})=h_{m}e^{-\frac{(\mathbf{x}-\mathbf{x}_{0})^{2}}{2r^{2}}}$, and
a horizontal scale, $r$, randomly chosen from an exponential probability
distribution, $f_{r}$, with a typical scale $r_{0}$ ,
$f_{r}(r)=\frac{1}{r_{0}}e^{-r/r_{0}}$, equivalent to the distribution of
circle radii in the void model of Popović et al. [14]. To prevent having
unrealistically narrow and high mounds, we prescribed the height of each
mound, $h_{m}$, to be proportional to its horizontal scale. Optionally, the
mounds may also be elongated along a certain axis and the axes may be
preferentially aligned to simulate anisotropy in the snow dunes. Including
anisotropy did not change any of the conclusions of our model. The “snow dune”
surface has a height distribution that is well-fit with a gamma distribution
(see Popović et al. [15] for more details). Parameters of the surface height
distribution depend on the density of mounds placed, $\rho\equiv
N_{\text{m}}\frac{r_{0}^{2}}{L^{2}}$, where $N_{\text{m}}$ is the number of
mounds placed within the domain of size $L$. When few mounds are placed, such
that $\rho\ll 1$, the height distribution is highly right-skewed, while when
many mounds are placed, such that $\rho\gg 1$, the height distribution
converges to a Gaussian. For this reason the percolation threshold also
depends on the density of mounds. Specifically, when $\rho$ is large,
$p_{c}=0.5$ due to the symmetric height distribution. When the height
distribution significantly deviates from a Gaussian, $p_{c}<0.5$. In Popović
et al. [15], we find that the “snow dune” topography with approximately
$0.2<\rho<0.5$ reproduces the LiDAR measurements of snow on first-year ice.
Using $\rho$ in this range, we find that the percolation threshold of a
realistic “snow dune” topography lies approximately between 0.4 and 0.5.

In addition to the “snow dune” topography, we also used two other types of
topographies. The first type we called the “diffusion” topography (Fig. S1b).
Ponds that form on this surface are shown in Fig. S1e. This topography is
generated by first independently assigning a height to each grid point at
random, and then letting this configuration diffuse for some time. Thus, to
generate this topography, we numerically solve the diffusion equation.
Diffusion smooths out variability in height up to a certain length scale that
depends on the time allowed for diffusion. This length scale then determines
the length scale $l_{0}$. To obtain the desired $l_{0}$, we set the diffusion
coefficient to 1 and tune the time for diffusion. This surface has a Gaussian
height distribution. Due to this fact, the percolation threshold is
$p_{c}=0.5$. We note that a similar type of topography has been previously
used to model melt pond geometry [2].

We called the other type of topography we used to test our model the
“Rayleigh” topography (Fig. S1c). Ponds that form on this surface are shown in
Fig. S1f. This topography is made by generating two “diffusion” topographies,
$h_{1}(\mathbf{x})$ and $h_{2}(\mathbf{x})$, with the same length-scale but
initialized with a different random configuration, and then taking the square
root of the sum of their squares,
$h(\mathbf{x})=\sqrt{h_{1}^{2}(\mathbf{x})+h_{2}^{2}(\mathbf{x})}$. This
surface has a non-symmetric Rayleigh height distribution (hence the name). The
percolation threshold on this surface is $p_{c}\approx 0.4$. A percolation
threshold that deviates from 0.5 is the main reason we tested this surface.
“Diffusion” and “Rayleigh” surfaces generated as above have no obvious
relation to any physical feature in sea ice, but are useful in order to test
the universality of drainage.

Figure S1: a-c: Examples of different types of topographies we used. Red
colors indicate highs of the topography, while blue colors indicate lows of
the topography. a) A “snow dune” topography. b) A “diffusion” topography. c) A
“Rayleigh” topography. d-f: Examples of ponds on different topographies. d)
Ponds on a “snow dune” topography. e) Ponds on a “diffusion” topography. f)
Ponds on a “Rayleigh” topography. g) A binarized image of a real melt pond
photograph taken on August 14th during the HOTRAX mission.

When generating each of the topographies above, for convenience, we assumed
periodic boundary conditions. To estimate $p_{c}$ for each synthetic
topography, we cut the topography with a horizontal plane and shifted it up or
down, tracking the connected clusters of the surface that lay below this
plane. We found a position of the plane at which a cluster that spans the
domain first appears, and estimated $p_{c}$ as the fraction of the surface
below the plane at this level. A constant $c$ in the parameter
$\eta=cNl_{0}^{2}/L^{2}$ is a property of a topography type in a similar way
as $p_{c}$. In section S6, we show that this constant can be related to other
properties of the surface such as the percolation threshold and the amplitude
of the cluster correlation length near the percolation threshold. However,
here we determined $c$ empirically by finding a value for which curves
$\Pi=g(\eta)$ on synthetic topographies best fit $g(\eta)$ estimated
analytically using Eq. LABEL:eqn:4 of the main text according to the least-
squares metric. Similar to $p_{c}$, a constant $c$ for “snow dune”
topographies could in principle depend on the density of mounds, $\rho$.
However, we found that $c\approx 3$ for all tested densities of mounds. We
found that for the “diffusion” surface $c\approx 4.1$ and $c\approx 3$ for the
“Rayleigh” surface.

##  S4 Comparing the full 2d model to the estimate

In Fig. LABEL:fig:3 of the main text, we compared the full 2d hole model to
Eqs. LABEL:eqn:12 to LABEL:eqn:14 of the main text. To make this comparison,
we needed to specify the hole opening distribution, $F$, the universal
function $g$, and the parameters $T_{\text{h}}$, $T_{\text{m}}$, $t_{0}$,
$l_{0}$, $L$, $N_{0}$, $c$ and $p_{c}$ in Eqs. LABEL:eqn:12 to LABEL:eqn:14 of
the main text. Here, we describe in detail how we related the parameters
inputed in the full 2d model to parameters that enter Eqs. LABEL:eqn:12 to
LABEL:eqn:14 of the main text. Note that the parameters of the 2d hole model
are not related to Eqs. LABEL:eqn:12 to LABEL:eqn:14 via Eqs. LABEL:eqn:5 to
LABEL:eqn:11 of the main text that describe how measurable physical parameters
enter the pond evolution equations since not all of these physical parameters
are prescribed in the 2d model (see Table S1 for the list of parameters
directly prescribed in the 2d model).

  1. 1.

We derived the hole opening distribution, $F$, as the cumulative distribution
of the distribution of critical temperatures, $f_{\theta}$, defined in the
full 2d model. This relationship is correct because we used a constant rate of
bulk ice temperature increase in the 2d model. Since, as a default, we used a
normal $f_{\theta}$, the default $F$ was a cumulative normal distribution.

The universal function, $g$, is simply the solution to Eq. LABEL:eqn:4 of the
main text.

  2. 2.

The hole opening timescale, $T_{\text{h}}$, was directly prescribed in the
full 2d model as the time for the bulk ice temperature to increase by one
standard deviation of the distribution $f_{\theta}$. Therefore, we used this
timescale in Eqs. LABEL:eqn:12 to LABEL:eqn:14 of the main text as well.

  3. 3.

The memorization timescale, $T_{\text{m}}$, was not prescribed directly in the
2d model, so we had to derive it in order to use it in Eqs. LABEL:eqn:12 to
LABEL:eqn:14 of the main text. To find $T_{\text{m}}$ in simulations with
$\text{d}H/\text{d}t=0$, we used ice thickness, $H$, and the preferential
ponded ice melt rate, $\text{d}h_{\text{diff}}/\text{d}t$, prescribed in the
full 2d model to simultaneously solve for $p_{\text{min}}$ and $T_{\text{m}}$
using Eqs. LABEL:eqn:7 and LABEL:eqn:13 of the main text.

  4. 4.

The center of the hole distribution, $t_{0}$, fluctuates slightly for each run
of the 2d model. For this reason in Eqs. LABEL:eqn:12 to LABEL:eqn:14 of the
main text, we did not use the approximate relation for $t_{0}$ defined in Eq.
LABEL:eqn:9 of the main text. Instead, we explicitly recorded the timing of
the center of the hole distribution relative to the timing of the first hole
for each run of the 2d model and used that number in the corresponding
estimate. Using Eq. LABEL:eqn:9 of the main text to find $t_{0}$ also gives
relatively good results, but there are noticeable deviations between the
estimate and the actual model solution.

  5. 5.

For the typical pond size, $l_{0}$, we used the length at which the
autocorrelation function for a configuration of ponds in the 2d model drops by
a factor of $e$.

  6. 6.

Domain size, $L$, is simply the number of grid points on the side of the
domain in the full 2d model.

  7. 7.

Potential number of holes, $N_{0}$, in simulations where each grid point
represents a potential hole was equal to the total number of grid points,
$N_{0}=L^{2}$. In simulations where there can be more than one hole per grid
point, $N_{0}$ is correspondingly larger.

  8. 8.

Parameters $c$ and $p_{c}$ are properties of the topography type, so in the
estimates, we used the parameters of the corresponding topographies we used in
the full 2d simulations. As we described in section S3, for each topography
type, we determined $c$ by least-squares fitting the drainage curves to Eq.
LABEL:eqn:4 of the main text, and we estimated $p_{c}$ as the fraction of the
surface below a level plane that cuts through the topography such that a
connected level set that spans the domain first appears.

Relating the parameters as described above, we were able to uniquely compare
estimates using Eqs. LABEL:eqn:12 and LABEL:eqn:13 of the main text to the
full 2d hole model for stage II in the absence of ice thinning. This
comparison was shown in Fig. LABEL:fig:3a of the main text. In Fig.
LABEL:fig:3b of the main text, we included ice thinning in the full 2d model
to simulate stage III. There, we prescribed a constant ice thinning rate,
$\text{d}H/\text{d}t$, and the initial ice thickness, $H_{0}$, such that ice
thickness evolves as

$H(t)=H_{0}-\frac{\text{d}H}{\text{d}t}t\text{ . }$ (S3)

When making the estimate with ice thinning included, most of the parameters
that enter Eqs. LABEL:eqn:12 to LABEL:eqn:14 of the main text can be related
to the 2d model parameters in the same way as described above, with the most
notable difference being the memorization timescale, $T_{\text{m}}$. Equation
LABEL:eqn:7 of the main text defines $T_{\text{m}}$, as the time for ponded
ice to melt below sea level in terms of ice thickness and pond coverage. Since
$T_{\text{m}}$ is the time it takes for ponded ice to melt through the
thickness of post-drainage freeboard, the thickness and the pond coverage that
enter Eq. LABEL:eqn:7 of the main text represent the post-drainage parameters.
Therefore, to get $T_{\text{m}}$, we need to solve Eqs. LABEL:eqn:7 and
LABEL:eqn:13 of the main text, and Eq. S3 simultaneously for post-drainage
coverage, $p_{\text{min}}$, memorization timescale $T_{\text{m}}$, and post-
drainage thickness, $H(T_{\text{m}})$.

In addition to predicting the time at which pond bottoms fall below sea level,
$T_{\text{m}}$ estimated as above also coincides fairly accurately with the
time at which the water level throughout the domain reaches sea level. This is
because, if pond bottoms lay below sea level, ponds do not become disconnected
as the drainage progresses, so a single hole can drain an entire pond, making
the drainage highly efficient after $T_{\text{m}}$. Note that this is only
true because, in our model, ice at the beginning of stage II is flat compared
to the freeboard thickness, so that pond bottoms are also flat and all of
ponded ice falls below sea level at roughly the same time. If this were not
the case, it could happen that drainage through a hole splits a pond into
disconnected parts thereby leaving a part undrained, even if the pond bottom
was initially below sea level. The fact that drainage through holes becomes
efficient once pond bottoms fall below sea level means that Eq. LABEL:eqn:12
of the main text is approximately valid up to $T_{\text{m}}$, and that Eq.
LABEL:eqn:14 of the main text is approximately valid after $T_{\text{m}}$.
$T_{\text{m}}$ estimated in this way is also a crossover time when $p(t)$
estimated using Eq. LABEL:eqn:12 of the main text becomes less than $p(t)$
estimated using Eq. LABEL:eqn:14 of the main text. Like Eq. LABEL:eqn:13, Eq.
LABEL:eqn:14 of the main text cannot be solved directly since there the pond
coverage $p(t)$ depends on $T_{\text{m}}(t)$, which in turn depends on $p(t)$
through Eq. LABEL:eqn:15 of the main text. Therefore, to solve for pond
coverage evolution during stage III, we have to simultaneously solve for
$p(t)$ and $T_{\text{m}}(t)$ using Eqs. LABEL:eqn:14 and LABEL:eqn:15 of the
main text where $H(t)$ is given by Eq. S3.

##  S5 Satellite and field data analysis

In the main text we used field measurements and satellite data to constrain
our model and to test its predictions. In this section we describe these
measurements and discuss how we used them to estimate parameters that enter
our model.

###  S5.1 Estimating parameters based on field measurements

Polashenski et al. [11] collected extensive field data aimed at understanding
the formation and evolution of melt ponds. The experiment was performed on
land-fast first-year ice near Barrow, Alaska and was repeated during the
summers of 2008, 2009, and 2010, with the most extensive studies done during
2009. Measurements of pond coverage shown in Fig. LABEL:fig:4a of the main
text were made along a 200 m-long transect every few days during 2009.
Additionally, they made measurements including ice thickness, temperature and
salinity, ice and pond albedo, meltwater production and drainage, and ice and
snow topography. Based on their measurements and data collected during other
studies, we were able to estimate most of the parameters that enter the pond
evolution equation. These values are reported as default values in Table
LABEL:tab:1 of the main text along with a plausible range when it was possible
to estimate that.

We estimated the linear size of the drainage basin, $L$, based on the remark
of Polashenski et al. [11] that the drainage basin where the measurements were
made is 1.3 km$\times$1.7 km in size. Therefore, we used $L\approx 1.5\text{
km}$. We estimated $T_{\text{m}}$ using Eq. LABEL:eqn:7 of the main text
coupled with their measurements of ice thickness ($H\approx 1.2\text{ m}$),
pond and ice albedo ($\alpha_{p}\approx 0.25$, $\alpha_{i}\approx 0.6$), and
measurements of solar flux from the nearby weather station
($F_{\text{sol}}\approx 254\text{ W}\text{ m}^{-2}$ average flux during field
experiment). We estimated the total number of potential holes, $N_{0}$, based
on measurements of Golden [6] who found that the density of brine channels was
between 60 and 120 per $\text{m}^{2}$ in their measurements. We thus used
$N_{0}\approx L^{2}100\text{ m}^{-2}$. We note that depending on the
conditions during ice growth, the density of brine channels may vary
substantially [20]. Nevertheless, as we show in section S8, pond coverage
depends only weakly on brine channel density, so even an order of magnitude
difference in brine channel density leads to only several percent change in
pond coverage.

Polashenski et al. [11] did not estimate the pond size $l_{0}$. However, based
on Popović et al. [15] the pond size seems to be constant between different
years with $l_{0}\approx 5.5\text{ m}$. Popović et al. [14] estimated the
percolation threshold for late-summer ponds for two different years and found
it to be around 0.3 and 0.4 for these years. Our simulations on the “snow
dune” topography predict $0.4<p_{c}<0.5$. Based on these estimates and
simulations we determine that the percolation threshold likely lays somewhere
between 0.3 and 0.5. The constant $c$ had only a small effect in our
simulations so we chose $c=3$, consistent with the snow dune topography. We
estimated the center of the hole distribution, $t_{0}$, based on Eq.
LABEL:eqn:9 of the main text and we found $t_{0}\approx 13.4\text{ days}$.

Finally, we estimated $T_{\text{h}}$ using Eqs. LABEL:eqn:10 and LABEL:eqn:11
of the main text. This estimate is approximate for two reasons - 1) we are
uncertain about the mechanism that drives hole opening and 2) even if the
mechanism is correct, some of the parameters that enter Eqs. LABEL:eqn:10 and
LABEL:eqn:11 of the main text are difficult to estimate. In addition to
measurements of thickness, albedo and solar flux, we used a measurement by
Polashenski et al. [11] of ice salinity, $S=3\text{ ppt}$. The extinction
coefficient, $\kappa$, is relatively well-documented and we used
$\kappa=1.5\text{ m}^{-1}$ [19]. Ice conductivity, $k$, can vary because ice
and brine have different conductivities. Using the same conductivity
parameterization as in Polashenski et al. [11], and assuming a brine volume
fraction between 0 and 0.5, we found that conductivity can vary between 1.3
and 2$\text{ W}\text{ m}^{-1\circ}\text{C}^{-1}$. We used $k=1.8\text{
W}\text{ m}^{-1\circ}\text{C}^{-1}$ which corresponds to a brine volume
fraction of 0.1. Estimating parameters $\theta_{0}$, $\Delta\theta$, $c^{*}$,
and $z^{*}$ is difficult. To estimate the reference temperature, $\theta_{0}$,
and the temperature range, $\Delta\theta$, we used observations made by
Polashenski et al. [12]. They noted that holes tend to begin to open when ice
interior temperature reaches roughly $-1.6^{\circ}$C and open completely when
ice temperature reaches around $-0.9^{\circ}$C. Therefore, we chose the
temperature range to be $\Delta\theta=0.7^{\circ}\text{C}$. Based on these
observations and our solutions to the full heat equation, we chose the
reference temperature to be $\theta_{0}=-1.2^{\circ}$C, and we used the shape
parameter, $c^{*}=2$ (see section S7). Polashenski et al. [12] found that
freshwater plugs do not form at a single depth within the ice, but rather
throughout a range of depths. Nevertheless, we used a single $z^{*}$ for
simplicity. Polashenski et al. [12] took photographs that show that freshwater
plugs tend not to form in the upper 0.6 m of the ice, and we estimated roughly
$z^{*}=0.6\text{ m}$ for this reason. Future models will need to include the
fact that tracking a single layer within the ice is insufficient to capture
the process of hole formation.

###  S5.2 Comparing the estimated pond coverage to Polashenski et al. [11]

With the physical parameters estimated above, we were able to find all of the
parameters that enter Eqs. LABEL:eqn:12 to LABEL:eqn:14 of the main text using
Eqs. LABEL:eqn:7 and LABEL:eqn:9-LABEL:eqn:11 of the main text. We used these
solutions to compare our estimate to pond coverage evolution observed by
Polashenski et al. [11] in Fig LABEL:fig:4a of the main text. Parameters
$\theta_{0}$, $\Delta\theta$, $c^{*}$, and $z^{*}$ clearly have larger
uncertainty than other parameters that can be directly measured. Therefore, to
find the wider error bars in Fig LABEL:fig:4a of the main text (light shaded
region), we assumed a $50\%$ error in the parameters $\theta_{0}$,
$\Delta\theta$, $c^{*}$, and $z^{*}$, and used a $10\%$ measurement error for
all other parameters. For narrower error bars in Fig LABEL:fig:4a of the main
text (dark shaded region) we used the $10\%$ measurement error for all
parameters. These measurement errors are not representative of real
uncertainties, and we chose them mainly for demonstrative purposes. In fact,
these errors are likely still underestimated. For this reason, our estimates
of pond coverage are uncertain and all we can say is that they are at least
consistent with our model.

###  S5.3 MODIS satellite data

A dataset of pond coverage fraction estimates based on satellite measurements
is available for download at https://cera-www.dkrz.de/WDCC/ui/cerasearch/. In
this dataset, pond fraction is estimated using an artificial neural network
that takes as input the reflectances from three channels of the visible
spectrum of the MODIS instrument aboard the Earth Observation Satellite TERRA
[16]. The observations are available every 8 days, have a spatial resolution
of 12.5km, and cover the entire Arctic. To make Fig. LABEL:fig:4b of the main
text, we collected estimated pond fractions for all grid-cells at all times
during a particular year for which there were measurements available. We only
considered measurements that showed a non-zero pond coverage. We then
summarized these data as probability distributions for several different
years. Due to the relatively coarse resolution of the measurements, the
estimated maximum pond coverage is likely somewhat underestimated.

##  S6 Deriving the universal drainage curve

Here, we will motivate a form of the universal function $g(\eta)$,
highlighting the origin of its universality. The goal is to estimate the
change in pond coverage fraction after a hole opens at a random location on
the ice surface. We are assuming that the first several holes have already
driven the ponds to the percolation threshold.

We introduce the cluster correlation function, $G(l)$, as the probability
that, given a randomly chosen point on a pond, another point of distance $l$
away will belong to the same pond. Equivalently, $G(l)$ can be viewed as the
average fraction of a circle of radius $l$ that belongs to the same pond as a
point chosen randomly on any ponded location. Therefore, a quantity $G(l)2\pi
l\text{d}l$ is an average area of a ring of radius $l$ and thickness
$\text{d}l$ connected to a hole opened randomly at some ponded location, and
an integral of $G$ over the entire 2d domain gives the average area of a pond
connected to such a hole. We need to multiply this integral by the pond
coverage fraction, $p$, to account for the fact that in our model holes can
open anywhere on a surface, not just within ponds. Therefore the mean fraction
of the surface connected to a randomly placed hole is

$\Delta p=p\frac{2\pi\int_{0}^{\infty}G(l)l\text{d}l}{L^{2}}\text{ . }$ (S4)

Equation S4 would represent the average change in coverage fraction after
opening a hole, if that hole drained the entire pond it was connected to.
However, a hole can only drain a fraction of an entire pond. Let us denote
that fraction as $f_{\text{drain}}$. The fraction of the pond drained depends
on the relative size of the pond to the autocorrelation length, $l_{0}$.
Namely, when a pond has a size comparable to $l_{0}$, a hole can drain a
significant fraction of that pond. On the other hand, when a pond is very
large compared to $l_{0}$, a small amount of drainage would quickly lead to
the formation of disconnected regions, and a hole would be unable to drain a
significant fraction of the pond. We can describe a typical linear size of the
largest ponds using a “cluster correlation length,” $\xi$, defined roughly as
the linear extent of the largest connected pond or, equivalently, as the
length beyond which $G(l)$ rapidly falls to 0. Therefore, a fraction
$f_{\text{drain}}$ can be viewed as a function of $l_{0}/\xi$ such that
$f_{\text{drain}}(l_{0}/\xi)\rightarrow 0$ when $l_{0}/\xi\rightarrow 0$ and
$f_{\text{drain}}(l_{0}/\xi)\rightarrow 1$ when $l_{0}/\xi\rightarrow\infty$.
The correlation length, $\xi$, depends on the deviation of pond coverage from
the percolation threshold since connected ponds become ever larger as the pond
coverage approaches $p_{c}$. Thus, close to the percolation threshold, ponds
are large and $l_{0}/\xi\ll 1$, so, based on a general Taylor expansion and
the above limiting behavior, we have

$f_{\text{drain}}(\frac{l_{0}}{\xi})\approx
a\frac{l_{0}}{\xi}+\text{...}\quad\text{ for
}\quad\big{(}1-\frac{p}{p_{c}}\big{)}\ll 1\text{ ,}$ (S5)

where $a$ is some constant. We note that this expansion hides an implicit
assumption that $f_{\text{drain}}(l_{0}/\xi)$ is analytic in the limit
$l_{0}/\xi\rightarrow 0$ and that $a$ is non-zero. In principle, a form
$f_{\text{drain}}(l_{0}/\xi)\propto(l_{0}/\xi)^{\alpha}$ as
$\l_{0}/\xi\rightarrow 0$ for any $\alpha>0$ would also be possible. So, we
can justify the above expansion only a posteriori, by showing that our
theoretical prediction matches the full 2d simulations well. Based on the
entire discussion above we can now estimate the pond fraction drained after
opening a hole as

$\frac{\text{d}p}{\text{d}N}=-p\frac{2\pi\int_{0}^{\infty}G(l)l\text{d}l}{L^{2}}f_{\text{drain}}(\frac{l_{0}}{\xi})\text{
. }$ (S6)

Each term on the right-hand side is a function of pond coverage, so this
equation defines a curve of pond coverage as a function of the number of open
holes.

Figure S2: a) Comparing $g(\eta)$ estimated by running a single 2d simulation
on a synthetic “snow dune” topography with no melt (black line) with a
solution of Eq. S6 using the full expressions Eqs. S11 to S13 (solid red line)
and a solution to universal Eq. S10 that assumes asymptotic forms of $G(l)$
and $\xi$, Eqs. S7 and S8 (dashed red line).

To close Eq. S6, we need to estimate the integral
$\int_{0}^{\infty}G(l)l\text{d}l$ and express $\xi$ in terms of $p$. General
percolation theory shows that in the limit of infinite domain size,
$L\rightarrow\infty$, close to the percolation threshold, $(1-p/p_{c})\ll 1$,
and at distances larger than the characteristic scale, $l\gg l_{0}$, $G(l)$
and $\xi$ behave as

$\displaystyle G(l)=p(\frac{l}{l_{0}})^{-5/24}e^{-l/\xi}\text{ , }$ (S7)
$\displaystyle\xi=Al_{0}(1-\frac{p}{p_{c}})^{-4/3}\text{ , }$ (S8)

for any model within the percolation universality class [7]. Here, $5/24$ and
$4/3$ are universal exponents that come from percolation theory, and $A$ is a
non-universal order one constant. At distances smaller than $l_{0}$, $G(l)$
depends on the shape of individual ponds, and, thus, on the details of the
topography, but is nevertheless constrained to go to 1 as $l\rightarrow 0$.
The facts that 1) both the full $G(l)$ and its asymptotic form, Eq. S7, have
finite integrals as $l\rightarrow 0$, and 2) that $G(l)$ is significantly
above 0 for $l<\xi$, mean that, if $\xi$ is large compared to $l_{0}$, the
contribution of short lengths, $l\ll l_{0}$, to the integral
$\int_{0}^{\infty}G(l)l\text{d}l$ can be ignored for both the full $G(l)$ and
its asymptotic form, so we can use the asymptotic form, Eq. S7, to perform the
integration. If the domain size is not infinite, the size of the largest
connected ponds is limited by $L$, and, therefore, $\xi$ will not diverge as
predicted by Eq. S8, but will instead be limited by the domain size. Based on
the above discussion, we can use $G(l)$ and $\xi$ given by Eqs. S7 and S8 to
close the pond drainage equation, Eq. S6, when $l_{0}\ll\xi\ll L$. Moreover,
to ensure that we can use Eqs. S7 and S8 with Eq. S6 we additionally have to
assume that the hole opening process does not disrupt the geometry of the
ponds, so that $G$ and $\xi$ have the same form as more and more holes open.
Again, the assumption of pond geometry remaining unchanged can only be
justified a posteriori, by direct comparison with full 2d simulations.

We now show that pond drainage is universal in the limit $l_{0}\ll\xi\ll L$,
i.e., when the pond coverage is neither too close nor too far from the
percolation threshold. In that case, we can use the asymptotic forms of $\xi$
and $G(l)$ and the asymptotic expansion Eq. S5 for $f_{\text{drain}}$. We can
explicitly perform the integral of $G$ over $l$ to get

$\frac{\text{d}p}{\text{d}N}=-C\frac{l_{0}^{2}}{L^{2}}p^{2}(1-\frac{p}{p_{c}})^{-19/18}\text{
, }$ (S9)

where $C=2\pi\Gamma(\frac{43}{24})aA^{19/24}$ is a numerical constant that
depends on $A$, $a$, and constants of integration. Defining
$\Pi\equiv\frac{p}{p_{c}}$ and $\eta\equiv c\frac{l_{0}^{2}}{L^{2}}N$, where
$c\equiv p_{c}C$, we recover Eq. LABEL:eqn:4 of the main text

$\frac{\text{d}\Pi}{\text{d}\eta}=-\Pi^{2}(1-\Pi)^{-19/18}\text{ . }$ (S10)

The appropriate boundary condition for this equation is that ponds are at the
percolation threshold when only a finite number of holes are open,
$\Pi(\eta=0)=1$, where we used the fact that $\eta=0$ if the domain is
infinite and $N$ is finite. We can see that Eq. S10 defines a universal
function $\Pi=g(\eta)$, since there are no non-universal factors that enter
it. From the above derivation, we see that the universality of $g(\eta)$ comes
from the universality of the percolation theory in the limit $l_{0}\ll\xi\ll
L$. This derivation also clarifies that the factor $c$ that enters $\eta$
depends on the type of topography because it is a combination of factors $a$,
$A$, and $p_{c}$ which all potentially depend on the details of the
topography. A solution to Eq. S10 is shown with a red dashed line in Fig. S2.
It agrees well with the simulations apart from at low pond coverage and at
pond coverage close to the percolation threshold. Discrepancy between this
solution and the simulations near $p_{c}$ is due to the fact that in this
region $\xi\sim L$, and the universality breaks down. We find that this region
of disagreement shrinks as we increase our domain size $L$.

We will now show that we can closely match the simulated $g(\eta)$ also in the
regime $p\rightarrow p_{c}$ if we use $G(l)$ that is valid for all $l$
(including $l<l_{0}$) and $\xi$ that takes into account the effect of domain
size. The matching in this case, however, will not be universal.

As we have discussed above, Eqs. S7 and S8 are asymptotic limits only valid at
large separations, $l$ (so that $l\gg l_{0}$), and when the pond coverage is
neither too close nor too far from $p_{c}$ (so that $l_{0}\ll\xi\ll L$). In
the Supplementary Information of Popović et al. [14], we have shown that the
$G(l)$ and $\xi$ that fit the measured pond statistics at all lengths, $l$,
and finite domain size, $L$, are

$\displaystyle
G(l)=(e^{-l/l_{0}}(1-p)+p)(1+\frac{l}{l_{0}})^{-5/24}e^{-l/\xi}\text{ , }$
(S11) $\displaystyle\xi=\xi_{\infty}(1-e^{-BL/\xi_{\infty}})\text{ , }$ (S12)
$\displaystyle\xi_{\infty}=Al_{0}(1-\frac{p}{p_{c}})^{-4/3}\text{ , }$ (S13)

where $\xi_{\infty}$ is the cluster correlation length if the domain size were
infinite and has the same form as Eq. S8, and $B$ is a non-universal order one
constant. These equations respect the limits for $l\gg l_{0}$ and
$L\rightarrow\infty$, and also ensure that $G(l\rightarrow 0)\rightarrow 1$
and that $\xi\propto L$ when $\xi_{\infty}\rightarrow\infty$. However, the
particular form of these functions is arbitrary, and we chose them due to
their simplicity and the fact that they match the measured $G(l)$ and $\xi$
highly accurately outside the scaling limit (see SI of Popović et al. [14]).
Even though these equations describe $G(l)$ and $\xi$ well for all $l$ and
$L$, they only hold close to the percolation threshold, i.e., when $\xi\gg
l_{0}$, so we still expect deviations from them when $p\ll p_{c}$ (i.e., when
$\xi\sim l_{0}$).

Finally, we show the solution to Eq. S6 using the full Eqs. S11 to S13. To be
able to evaluate $\text{d}p/\text{d}N$ in this case, we also need to assume a
concrete $f_{\text{drain}}$ consistent with the asymptotic behavior given in
Eq. S5 and the fact that we require $f_{\text{drain}}(l_{0}/\xi)\rightarrow 1$
when $l_{0}/\xi\rightarrow\infty$. To this end, we choose
$f_{\text{drain}}=1-e^{-l_{0}/\xi}$. This arbitrary choice works well because
of universality in the scaling limit that we discussed above. With this form
of $f_{\text{drain}}$, we can use Eqs. S11 to S13 to solve Eq. S6 and get a
function $p=p(N)$. This solution in coordinates $p/p_{c}$ and
$Nl_{0}^{2}/L^{2}$ is shown as a solid red line in Fig. S2. We can see that
the agreement with the curve arising from a 2d simulation is excellent apart
from a small discrepancy at small pond coverage. This discrepancy is due to
the breakdown of Eqs. S11 to S13 far from $p_{c}$, due to our ad-hoc choice of
$f_{\text{drain}}$, and due to the fact that opening many holes likely does
not leave pond geometry completely unchanged. However, even though we were
able to obtain a good match with simulations, this solution is not universal
under the rescaling $p\rightarrow\Pi$ and $N\rightarrow\eta$. Namely, changing
the non-universal constants $A$ and $B$ or $f_{\text{drain}}$ changes the form
of the solution in $\Pi$ and $\eta$ coordinates, so, to get a good match, we
had to tune $A$ and $B$ specifically. This again highlights the fact that the
universality occurs in the limit $l_{0}\ll\xi\ll L$, i.e., when ponds are much
larger then the characteristic scale of the topography (set e.g. by the size
of snow dunes), and still smaller than the size of the domain (set e.g. by the
ice floe size).

##  S7 The heat equation

In Eqs. LABEL:eqn:10 and LABEL:eqn:11 of the main text, we noted that the
order of magnitude of the hole opening timescale, $T_{\text{h}}$, can be
estimated based on physical properties of ice if we assume that internal
temperature at some depth within the ice is the relevant parameter that
determines when holes start to open. In this section, following Bitz and
Lipscomb [1], we will briefly explain how we arrived at this estimate.

Energy supplied to a unit mass of salty ice partly goes into warming the ice
and partly into melting the ice around brine pockets to dilute the brine and
bring it into equilibrium with the surrounding ice. Ono [10] showed that the
heat capacity of a unit mass of ice of salinity $S$ and temperature $\theta$
in degrees Celsius can be expressed as

$c(\theta,S)=c_{0}+\frac{\gamma S}{\theta^{2}}\text{ , }$ (S14)

where $c_{0}=2.11\text{ kJ}\text{ kg}^{-1\circ}\text{C}^{-1}$ is the heat
capacity of fresh ice and $\gamma=18\text{ kJ}\text{
kg}^{-1}\text{ppt}^{-1\circ}\text{C}$ is a constant. Equation S14 neglects the
typically small contribution from the fact that water in brine pockets has a
different heat capacity than ice. This equation shows that the heat capacity
increases with temperature and diverges when the temperature approaches the
melting point of fresh ice. During the summer, ice is in the process of
melting, so the ice interior is only several degrees below zero. With typical
salinities of 3 ppt or 4 ppt for first year ice, the salinity term in the heat
capacity during summer is roughly an order of magnitude larger than $c_{0}$,
so $c(\theta,S)\approx\frac{\gamma S}{\theta^{2}}$.

The ice interior warms because of heat diffusion and because sunlight
penetrates beneath the ice surface. This can be expressed with a partial
differential equation that determines the warming rate at depth $z$ within the
ice

$\rho_{i}c(\theta,S)\frac{\partial\theta}{\partial t}=\frac{\partial}{\partial
z}k\frac{\partial\theta}{\partial z}+\kappa F_{\text{0}}e^{-\kappa z}\text{ ,
}$ (S15)

where $\rho_{i}$ is the density of ice, $k$ is the thermal conductivity of the
ice, $\kappa$ is the extinction coefficient from Beer’s law, and
$F_{\text{0}}$ is the radiative flux that penetrates the upper surface of the
ice. The first term on the right hand side represents the contribution from
heat diffusion while the second term is the contribution from direct solar
heating at a depth $z$. Thermal conductivity, $k$, may in principle depend on
depth because brine and ice have different thermal conductivities. This
dependence is, however, typically weak. Equation S15 should be supplied with a
boundary condition that ice is at its salinity-dependent melting temperature
at the top and bottom surface. Furthermore, a temperature at the initial time
should also be specified. We solve this equation in Fig. S3 assuming a uniform
salinity profile throughout the ice column and ice bottom temperature fixed at
-1.7${}^{\circ}\text{C}$.

Figure S3:  a) Black lines represent represent the vertical temperature
profiles within the ice at different times obtained by solving Eq. S15. Curves
to the left correspond to earlier times. Colored horizontal lines represent
different choices of $z^{*}$ that are shown in panels b and c. b) Time
evolution of ice temperature at a fixed depth. Different colors stand for
different depths, $z^{*}$, shown in panel a. Horizontal dashed lines represent
the temperature at which the holes tend to start opening,
$\theta_{\text{min}}$, and the temperature at which a significant fraction of
brine channels are open, $\theta_{\text{max}}$. The temperature range,
$\Delta\theta$, is also marked. Vertical dashed lines represent times at which
the temperature at depth $z^{*}=0.6\text{ m}$ crosses $\theta_{\text{min}}$
and $\theta_{\text{max}}$, which defines the hole opening timescale,
$T_{\text{h}}$. A reference temperature, $\theta_{0}$, estimated as the middle
of the range between $\theta_{\text{min}}$ and $\theta_{\text{max}}$ is also
marked with a horizontal dashed line. c) Dependence of the warming rate on
current ice temperature. Different colors stand for different depths, $z^{*}$,
marked in panel a. Black dashed curves show estimates using Eq. S16 with
$c^{*}=2$ and treating $\theta_{0}$ as a variable. Horizontal solid colored
lines show $\text{d}\theta/\text{d}t$ estimated as
$\Delta\theta/T_{\text{h}}$, where $T_{\text{h}}$ is found numerically as in
panel b. Horizontal colored dotted lines show $\text{d}\theta/\text{d}t$
estimated using Eq. S16 with $c^{*}=2$ and a reference temperature
$\theta_{0}=0.5(\theta_{\text{max}}+\theta_{\text{min}})$.

As we have suggested in the main text, our strategy is to estimate
$T_{\text{h}}$ by estimating the warming rate
$\frac{\text{d}\theta}{\text{d}t}$ at some depth, $z^{*}$, within the ice
where ice plugs tend to form. In general, the warming rate in Eq. S15 cannot
be simply characterized as it depends on the full vertical profiles of
temperature and salinity and their history. Nevertheless, we can approximate
it as

$\frac{\text{d}\theta}{\text{d}t}\approx\frac{\theta_{0}^{2}}{\rho_{i}\gamma
S}\Big{(}c^{*}k\frac{\theta_{0}}{H^{2}}+(1-\alpha_{p})F_{\text{sol}}\kappa
e^{-\kappa z^{*}}\Big{)}\text{ . }$ (S16)

Here, we have assumed that the fresh ice heat capacity, $c_{0}$, is
negligible, so that $c(\theta,S)\approx\frac{\gamma S}{\theta^{2}}$. We have
also assumed that thermal conductivity, $k$, does not depend on brine volume
fraction, and we have estimated the amount of radiation that penetrates the
upper surface as $F_{\text{0}}\approx(1-\alpha_{p})F_{\text{sol}}$, where
$\alpha_{p}$ is the pond albedo and $F_{\text{sol}}$ is the solar radiative
flux. We only consider pond albedo because we are only interested in holes
that open beneath ponds. A constant $c^{*}$ accounts for the shape of the
temperature profile, while $\theta_{0}$ and $z^{*}$ are the reference
temperature and depth at which we are estimating the warming rate.

To be able to estimate the hole opening timescale, $T_{\text{h}}$, using Eq.
S16, we need to relate the parameters of that equation to physically
meaningful quantities. We can define $T_{\text{h}}$ as the time it takes for
ice at a particular depth to warm from a temperature, $\theta_{\text{min}}$,
at which only a small fraction of the brine channels can become holes, to a
temperature, $\theta_{\text{max}}$, at which a significant fraction of the
brine channels become holes (see Fig. S3b). We denote the difference
$\theta_{\text{max}}-\theta_{\text{min}}$ as $\Delta\theta$. We thus consider
$T_{\text{h}}$, $\theta_{\text{min}}$, $\theta_{\text{max}}$, and
$\Delta\theta$ to be measurable, physically meaningful quantities. We can
approximate $\theta_{0}$ in Eq. S16 to be the middle of the range between
$\theta_{\text{max}}$ and $\theta_{\text{min}}$,
$\theta_{0}\approx\frac{\theta_{\text{min}}+\theta_{\text{max}}}{2}$. We
consider $z^{*}$ to have physical significance as the depth at which ice plugs
form, and can thus be estimated using measurements. Finally, the parameter
$c^{*}$ takes into account the shape of the temperature profile within the ice
column. As such, it depends on a multitude of factors such as the boundary
conditions, the salinity profile, or the depth at which we are estimating
$\frac{\text{d}\theta}{\text{d}t}$. To keep matters simple, we choose
$c^{*}=2$, corresponding to a quadratic temperature profile with the top and
bottom ice temperatures fixed at 0${}^{\circ}\text{C}$.

To determine how reliable the approximations we made above are, we solved the
full heat equation, Eq. S15, varying the equation parameters. In Fig. S3c, we
show $\frac{\text{d}\theta}{\text{d}t}(\theta)$ found by solving Eq. S15
numerically, and compare it to $\frac{\text{d}\theta}{\text{d}t}(\theta_{0})$
estimated using Eq. S16 treating $\theta_{0}$ as a variable (compare colored
and dashed lines in Fig. S3c). We can see that in all cases, the estimate and
the numerical simulations behave qualitatively in the same way and are of
similar magnitude (within around 20% of each other). In Fig. S3c, we also show
$\frac{\Delta\theta}{T_{\text{h}}}$, where we estimate $T_{\text{h}}$ by
solving Eq. S15 numerically and compare it to
$\frac{\text{d}\theta}{\text{d}t}$ estimated using Eq. S16 with
$\theta_{0}=\frac{\theta_{\text{min}}+\theta_{\text{max}}}{2}$. Again, we find
that our approximations are of the same order of magnitude as the simulations
(in this case within 10% of each other), which implies that $T_{\text{h}}$ can
be reasonably estimated using Eq. S16.

For each particular numerical solution of Eq. S15, we can improve the match
between the simulations and estimates by tuning $c^{*}$ or changing the
definition of $\theta_{0}$. Depending on the depth, boundary conditions, shape
of the salinity profile, and the strength of direct solar heating, we find
that $c^{*}$ can be anything between roughly 1 and 10, while $\theta_{0}$ can
be modified to $\theta_{0}=(1-r)\theta_{\text{max}}+r\theta_{\text{min}}$,
with the numerical factor $r$ anywhere between 0.4 and 0.5. Nevertheless,
using $c^{*}=2$ and
$\theta_{0}=\frac{\theta_{\text{min}}+\theta_{\text{max}}}{2}$, as we did
above, and varying the parameters of the numerical simulation, we find that
$\frac{\text{d}\theta}{\text{d}t}(\theta_{0})$ estimated using Eq. S16 and
treating $\theta_{0}$ as variable is always of the same order of magnitude as
$\frac{\text{d}\theta}{\text{d}t}(\theta)$ estimated using Eq. S15 for all
$\theta_{0}$ and $\theta$, while $\frac{\text{d}\theta}{\text{d}t}$ estimated
at $\theta_{0}=\frac{\theta_{\text{min}}+\theta_{\text{max}}}{2}$ using Eq.
S16 is always of the same order of magnitude as
$\frac{\Delta\theta}{T_{\text{h}}}$, with $T_{\text{h}}$ estimated by solving
Eq. S15 numerically. This means that Eq. S16 with $c^{*}=2$ and $r=0.5$ is a
reasonable first-order approximation to the full heat equation under any
configuration, with errors in estimated $\frac{\text{d}\theta}{\text{d}t}$
likely being on the order of 10-20% of the full solution based on the results
in Fig. S3c.

Finally, using Eq. S16, we can estimate the hole opening timescale,
$T_{\text{h}}$ as

$T_{\text{h}}=\frac{\Delta\theta}{\text{d}\theta/\text{d}t}\text{ , }$ (S17)

with $\frac{\text{d}\theta}{\text{d}t}$ estimated using Eq. S16. In addition
to the fact that this estimate may disagree with the full solution to the heat
equation as we have discussed above, there are several additional reasons why
it is uncertain - 1) it is difficult to precisely define and determine
parameters $\theta_{\text{max}}$, $\theta_{\text{min}}$, and $z^{*}$ from
measurements, 2) it is unclear whether these parameters themselves depend on
physical properties such as ice salinity or thickness, and, most importantly,
3) it is unclear whether the ice interior temperature is in fact a good
indicator of when holes start to open. In particular, the assumption that ice
interior temperature is a good indicator of hole opening leads to some
counter-intuitive conclusions that we discuss in section S8 and section
LABEL:sec:discussion of the main text. Namely, it predicts that holes open
more slowly when ice is more saline or when it is closer to its melting point,
since in these cases $\frac{\text{d}\theta}{\text{d}t}$ is lower due to a
higher higher heat capacity stemming from a higher brine volume fraction. This
is counter-intuitive because it implies that it is more difficult to create
holes when ice is more porous and when more heat is expended into melting the
ice pores. For these reasons, improving our understanding of hole formation
physics is crucial for understanding melt pond evolution.

##  S8 Dependence of pond coverage on physical parameters

In this section, we will discuss how pond coverage depends on the measurable
properties of the ice in our model. In particular, we will focus on the pond
coverage minimum, $p_{\text{min}}$, found using Eqs. LABEL:eqn:7,
LABEL:eqn:9-LABEL:eqn:11, and LABEL:eqn:13 of the main text. In the entire
discussion below, we will use the cumulative normal hole opening distribution,
$F$, and the set of default parameters defined in Table LABEL:tab:1 of the
main text.

Figure S4: Post-drainage pond coverage, $p_{\text{min}}$, as a function of
physical parameters using Eqs. LABEL:eqn:7, LABEL:eqn:9-LABEL:eqn:11, and
LABEL:eqn:13 of the main text. In each panel, we are varying one parameter and
assuming all other parameters are the defaults specified in Table LABEL:tab:1
of the main text. We are assuming a cumulative normal hole opening
distribution, $F$. Different colors in all panels except e and l stand for
different values of the depth at which ice plugs tend to form, $z^{*}$. a)
Pond coverage as a function of the basin length-scale, $L$. Here, we are
keeping the density of brine channels, $n_{0}=N_{0}/L^{2}$, constant. b) Pond
coverage as a function of the typical pond length-scale, $l_{0}$. c) Pond
coverage as a function of the percolation threshold, $p_{c}$. d) Pond coverage
as a function of the brine channel density, $n_{0}=N_{0}/L^{2}$. e) Pond
coverage as a function of ice thickness, $H$. The red curve assumes that ice
plugs form at a constant depth $z^{*}=0.6$ m. The black dashed curve assumes
that ice plugs always form in the middle of the ice. f) Pond coverage as a
function of ice salinity, $S$. g) Pond coverage as a function of solar flux,
$F_{\text{sol}}$. h) Pond coverage as a function of pond albedo, $\alpha_{p}$.
i) Pond coverage as a function of bare ice albedo, $\alpha_{i}$. j) Pond
coverage as a function of the reference temperature, $\theta_{0}$. k) Pond
coverage as a function of the temperature range, $\Delta\theta$. l) Pond
coverage as a function of the depth where ice plugs form, $z^{*}$.

As a summary of our results in the main text, here we recapitulate how the
minimum pond coverage, $p_{\text{min}}$, relates to physical parameters. After
some manipulation of Eqs. LABEL:eqn:5 to LABEL:eqn:13 of the main text, we
find

$\displaystyle
p_{\text{min}}=p_{c}g\big{(}\eta_{0}F(\tau_{\text{m}}-\tau_{0})\big{)}\text{ ,
}$ (S18) $\displaystyle\eta_{0}\equiv cn_{0}l_{0}^{2}\quad\text{ ,
}\quad\tau_{\text{m}}\equiv\frac{T_{\text{m}}}{T_{\text{h}}}\quad\text{ ,
}\quad\tau_{0}\equiv-F^{-1}\big{(}\frac{1}{n_{0}L^{2}}\big{)}$ (S19)
$\displaystyle\tau_{\text{m}}\approx\underbrace{\frac{l_{m}}{(\alpha_{i}-\alpha_{p})F_{\text{sol}}}\frac{\rho_{w}-\rho_{i}}{\rho_{w}}\frac{H}{1-p_{\text{min}}}}_{T_{\text{m}}}\quad\underbrace{\frac{\theta_{0}^{2}\Big{(}c^{*}k\frac{\theta_{0}}{H^{2}}+(1-\alpha_{p})F_{\text{sol}}\kappa
e^{-\kappa z^{*}}\Big{)}}{\Delta\theta\rho_{i}\gamma S}}_{1/T_{\text{h}}}$
(S20)

where $\alpha_{i}$ is the ice albedo, $n_{0}\equiv\frac{N_{0}^{2}}{L^{2}}$ is
the density of brine channels, and we have introduced the non-dimensional
parameters $\eta_{0}$, the number of brine channels per characteristic area of
the surface, $\tau_{\text{m}}$, the time for ponds to melt through the
freeboard relative to the hole opening timescale, $T_{\text{h}}$, and
$\tau_{0}$, the time between the first hole and the center of the hole
distribution, $F$, relative to $T_{\text{h}}$. Equation S20, relates the non-
dimensional timescale $\tau_{\text{m}}$ to measurable parameters using Eqs.
LABEL:eqn:5, LABEL:eqn:10 and LABEL:eqn:11 of the main text. Note that
$\tau_{\text{m}}$ depends on $p_{\text{min}}$ through the “memorization”
timescale, $T_{\text{m}}$, so Eqs. S18 and S20 have to be solved
simultaneously for $p_{\text{min}}$ and $\tau_{\text{m}}$.

The minimum pond coverage, $p_{\text{min}}$, is small when the differential
melt is weak or when the drainage is efficient. A larger $\tau_{\text{m}}$
means that ponds melt quickly relative to the timescale of hole opening, so
$p_{\text{min}}$ decreases with $\tau_{\text{m}}$. In particular, in terms of
$T_{\text{m}}$ and $T_{\text{h}}$, when $T_{\text{m}}$ is large, melt is weak,
so $p_{\text{min}}$ decreases with $T_{\text{m}}$, and when $T_{\text{h}}$ is
large, holes open slowly, so $p_{\text{min}}$ increases with $T_{\text{h}}$.
More brine channels per characteristic area increase the drainage efficiency,
so $p_{\text{min}}$ decreases with $\eta_{0}$. A wider time separation between
the initial hole and the center of the hole opening distribution means that
ponds will spend more time in the tail of the distribution with few holes
opening, decreasing the drainage efficiency, so $p_{\text{min}}$ increases
with $\tau_{0}$. We now discuss the dependence of $p_{\text{min}}$ on each of
the physical parameters that enter $\eta_{0}$, $\tau_{\text{m}}$, and
$\tau_{0}$. We show these relationships in Fig. S4, and we discuss each of the
panels in the list below.

  1. 1.

In Fig. S4a, we show that post-drainage pond coverage increases with size of
the drainage basin, $L$. This is because a larger basin makes it more likely
for the first hole to open somewhere within the domain, initiating the
drainage stage earlier, and increasing $\tau_{0}$.

  2. 2.

In Fig. S4b, we show that minimum pond coverage decreases with the typical
pond size, $l_{0}$. Physically, this is because larger ponds make it easier to
drain the surface. Recall that it takes about one hole per pond of size
$l_{0}$ to significantly drain the surface. Therefore, complete drainage can
be achieved with fewer open holes if the individual ponds are larger.

  3. 3.

Figure S4c shows that minimum pond coverage increases with the percolation
threshold, $p_{c}$. Pond coverage throughout its evolution from the end of the
initial phase above the percolation threshold is approximately proportional to
$p_{c}$. It is not exactly proportional since $\tau_{\text{m}}$ also depends
on $p_{c}$, making the dependence of $p_{\text{min}}$ on $p_{c}$ slightly non-
linear.

  4. 4.

Figure S4d shows the dependence of $p_{\text{min}}$ on the brine channel
density, $n_{0}$. We can see that this dependence is relatively weak - the
range of $n_{0}$ measured by Golden [6], from about $60\text{ m}^{-2}$ to
about $120\text{ m}^{-2}$, leads to only around 1% change in pond coverage.
Even an order of magnitude change in $n_{0}$ would lead to less than a 10%
change in pond coverage. Increasing $n_{0}$ leads to two competing effects. On
the one hand, a higher $n_{0}$ increases $\eta_{0}$ by increasing number of
potential holes, while on the other hand, it also increases $\tau_{0}$ by
increasing the time before a significant fraction of holes open. These two
effects largely cancel each other. As we show in section S9, if the hole
opening distribution, $F$, has an exponential tail, $F(x)\propto e^{-|x|}$,
the two effects exactly cancel. Since in our case $F(x)\propto e^{-x^{2}}$ in
the tail, some dependence on $n_{0}$ is retained, but the dependence remains
weak.

  5. 5.

In Fig. S4e, we show the dependence of $p_{\text{min}}$ on ice thickness, $H$.
Ice thickness enters the non-dimensional parameter $\tau_{\text{m}}$ by
increasing both the memorization and the hole opening timescale,
$T_{\text{m}}$ and $T_{\text{h}}$. In particular, larger $H$ means that, on
the one hand, ponds need to carve deeper depressions, thereby increasing
$T_{\text{m}}$, while on the other hand, ice also warms more slowly, thereby
increasing $T_{\text{h}}$. Because $H$ affects $T_{\text{h}}$, the effect of
ice thickness depends on hole formation physics. In Fig. S4e we show two
scenarios. The red line shows a case where ice plugs form at a constant depth,
$z^{*}=0.6\text{ m}$. We can see that in this case pond coverage peaks at an
ice thickness around 0.5 m after which it decreases to 0 with increasing
thickness. The black dashed line shows a case where ice plugs always form in
the middle of the ice. The dependence on ice thickness in this case is
qualitatively different than the fixed depth scenario. After peaking at
$H\approx 0.25\text{ m}$, pond coverage drops to its minimum at $H\approx
1.5\text{ m}$, after which it keeps increasing with increasing $H$. The
complicated manner in which pond coverage depends on $H$ in both of these
scenarios reflects two competing effects. So, depending on whether the effect
on $T_{\text{m}}$ or $T_{\text{h}}$ is stronger, pond coverage either
increases (if the effect on $T_{\text{h}}$ dominates) or decreases (if the
effect on $T_{\text{m}}$ dominates) with $H$. We note that Skyllingstad et al.
[18] found that there was no clear relationship between ice thickness and
post-drainage pond coverage in their detailed 3-dimensional model of ice and
ponds.

  6. 6.

In Fig. S4f, we show that in our model $p_{\text{min}}$ increases with ice
salinity, $S$. In our model, this is because saltier ice warms more slowly
since a large fraction of energy is deposited in melting ice around brine
pockets rather than going to warming the ice. This dependence is somewhat
counterintuitive since a higher salt content means a larger brine volume
fraction and therefore a larger potential for hole formation. For this reason
it may be an artifact of our assumption that holes open when a certain
temperature is reached. We note that Skyllingstad et al. [18] did not find a
significant relationship between salt content and minimum pond fraction in
their model.

  7. 7.

Figure S4g shows that $p_{\text{min}}$ depends only weakly on solar flux,
$F_{\text{sol}}$. This weak dependence is explained by the two competing
effects that nearly cancel each other. Namely, more solar radiation allows
ponds to preferentially melt ice faster, decreasing $T_{\text{m}}$, while also
warming the ice interior faster, thereby decreasing $T_{\text{h}}$.

  8. 8.

In Fig. S4h, we show how post-drainage pond coverage depends on pond albedo,
$\alpha_{p}$. More reflective ponds prevent ice beneath ponds from warming
quickly, increasing $T_{\text{h}}$, and also lead to a smaller contrast
between bare ice and ponds, $\Delta\alpha$, increasing $T_{\text{m}}$ as well.
However, the effect on $T_{\text{m}}$ dominates and pond coverage decreases
with increasing $\alpha_{p}$.

  9. 9.

In Fig. S4i, we show how post-drainage pond coverage depends on ice albedo,
$\alpha_{i}$. More reflective ice leads to a larger contrast between bare ice
and ponds, decreasing $T_{\text{m}}$, which leads to a higher pond coverage.

  10. 10.

In Fig. S4j, we show how $p_{\text{min}}$ depends on the reference
temperature, $\theta_{0}$. This reference temperature is approximately the
temperature at which holes tend to become open, and it is difficult to define
more precisely without further investigation into the hole opening process. We
can see that increasing $\theta_{0}$ leads to a rapid increase in pond
coverage in our model. As we have discussed in section S7, the reference
temperature affects the hole opening timescale by changing the warming rate.
Namely, close to the melting point of the ice, heat capacity increases,
rapidly decreasing the warming rate and increasing $T_{\text{h}}$ (see Fig.
S3c). Apart from increasing heat capacity, the warming rate also decreases
since heat diffusion weakens when ice is warmer. As with salinity, this may be
seen as counterintuitive since warmer ice contains a high brine volume
fraction which would be expected to aid the formation of holes.

  11. 11.

In Fig. S4k, we show that $p_{\text{min}}$ increases with the temperature
range, $\Delta\theta$. This is due to the fact that it takes a longer time to
warm a larger amount, increasing $T_{\text{h}}$.

  12. 12.

In Fig. S4l, we show that $p_{\text{min}}$ increases with the depth at which
ice plugs form, $z^{*}$. As less sunlight can penetrate deeper within the ice,
ice warms slower, increasing $T_{\text{h}}$. A caveat here is that we have
assumed that the rate of heat diffusion does not change with depth. As heat
diffusion is likely strong near the bottom of the ice, this relation likely
holds so long as the depth is not close to the ice bottom.

We believe that the assumptions we have made lead to dependencies on $L$,
$l_{0}$, $p_{c}$, $n_{0}$, $F_{\text{sol}}$, $\alpha_{p}$, $\alpha_{i}$,
$\Delta\theta$, and $z^{*}$ that are qualitatively correct. Despite this, it
is difficult to make quantitative claims since some of these parameters, such
as $\Delta\theta$ and $z^{*}$, are hard to constrain. In most panels of Fig.
S4 we showed that pond coverage varies appreciably by changing the under-
constrained parameter $z^{*}$, highlighting the quantitative uncertainty. We
note that prior to this investigation, it was not recognized that geometric
parameters $L$, $l_{0}$, and $p_{c}$ can have any effect on pond coverage. The
dependence on parameters $H$, $S$, and $\theta_{0}$ seems problematic. Namely,
as we have shown, the dependence on $H$ is highly sensitive to hole formation
physics. The dependence on $S$ and $\theta_{0}$ is somewhat counter-intuitive,
indicating that the assumption that temperature determines the onset of hole
formation needs to be examined in more detail. The above analysis shows that,
as we have already stressed several times, we need to improve our
understanding of hole formation physics in order to understand how pond
coverage depends on physical parameters.

Finally, we note that it is somewhat puzzling that the post-drainage coverage
fraction, $p_{\text{min}}$, is neither 0 nor $p_{c}$ in a typical realistic
situation. Namely, the function $g(\eta)$ behaves as $g(\eta)\rightarrow 0$
for $\eta\gg 1$ and $g(\eta)\rightarrow 1$ for $\eta\ll 1$, so there is only a
limited range of $\eta$ for which $p_{\text{min}}\neq 0$ and
$p_{\text{min}}\neq p_{c}$. The parameter $\eta$ falls within this range only
when $T_{\text{m}}\sim T_{\text{h}}$, so a question then arises of why the
time to melt through the thickness of the ice freeboard and the time for a
significant number of holes to open are comparable. Partly, it must be because
both phenomena are related to ice melt, but nevertheless, the two timescales
fundamentally depend on different properties of the ice. Moreover, since there
are many physical parameters that determine the value of $\eta$, its value
should wildly fluctuate from one set of environmental parameters to another,
making the outcome $\eta\ll 1$ or $\eta\gg 1$ likely. Therefore, the fact that
we typically observe $p_{\text{min}}\neq 0$ and $p_{\text{min}}\neq p_{c}$
perhaps indicates that some of the physical parameters that control pond
evolution are in fact not independent as we have assumed.

##  S9 The hole opening distribution

In the main text, we introduced a hole opening distribution,
$F(\frac{t-t_{0}}{T_{\text{h}}})$, as a cumulative of some arbitrary
probability density function, $f(\frac{t-t_{0}}{T_{\text{h}}})$. Here, we will
describe some ways that this distribution affects pond evolution.

First, we discuss how a time distribution, $f$, can be related to a measurable
distribution of brine channel properties, $f_{\theta}$. In general, we may
assume that there exists some underlying probability distribution,
$f_{\theta}(\frac{\theta-\theta_{0}}{\Delta\theta})$, that describes the
fraction of the holes that open when the ice interior temperature (or some
other bulk ice property) increases from $\theta$ to $\theta+\text{d}\theta$.
Such a distribution could in principle be measured in the field or modeled in
some way. It is related to the time distribution, $f$, as

$f(\frac{t-t_{0}}{T_{\text{h}}})=f_{\theta}(\frac{\theta-\theta_{0}}{\Delta\theta})\frac{\text{d}\theta}{\text{d}t}\text{
. }$ (S21)

This equation is valid even when the ice warming rate is not constant and
provides a basis to relate the hole opening distribution that enters Eqs.
LABEL:eqn:12 and LABEL:eqn:14 of the main text to measurable sea ice
properties.

Next, we discuss the variability of the distribution center, $t_{0}$. In the
main text, we remarked that the center of the hole opening distribution,
$t_{0}$, is approximately (see Eq. LABEL:eqn:9 of the main text)

$t_{0}\approx-T_{\text{h}}F^{-1}\big{(}\frac{1}{N_{0}}\big{)}\text{ . }$ (S22)

We argued that, because the timing of the opening of the first hole
fluctuates, the time until a given fraction of holes, $F(0)$, opens also
fluctuates, and, consequently, so does $t_{0}$. Here we will discuss the
probability distribution of $t_{0}$ over an ensemble of runs. Let us assume
that the first hole opens at some temperature $\theta_{i}$. $t_{0}$ is then
the time it takes for ice to warm from $\theta_{i}$ to $\theta_{0}$. There are
$N_{0}$ holes and each one opens at a different critical temperature
independently drawn from a distribution $f_{\theta}$. Thus, $\theta_{i}$ is a
minimum of a set of $N_{0}$ independent random variables drawn from a
distribution $f_{\theta}$. Therefore, the probability of finding some
$\theta_{i}$, and a corresponding $t_{0}$, is governed by extreme value
statistics (see, e.g., Coles [3] for an introduction to extreme value
statistics).

In general, depending on the tail of an underlying distribution, we expect the
probability of finding some minimum value from a large sample of independent
random variables to fall into one of three universal distributions. In our
case, we expect that $f_{\theta}$ has a relatively well-defined width. A
distribution with a well-defined width will likely have an exponential tail,
and such distributions have extreme value statistics described by a Gumbel
distribution

$f_{\text{Gumbel}}(x)=\frac{1}{\beta}e^{-(\frac{x-\mu}{\beta}+e^{-\frac{x-\mu}{\beta}})}\text{
, }$ (S23)

where $\mu$ is the mode of the distribution and $\beta$ is the scale.
Therefore, we can expect that, for each run, $t_{0}$ is drawn from this
distribution with a mode that is close to our approximation in Eq. S22.
Exactly how parameters $\mu$ and $\beta$ depend on $N_{0}$ and $T_{\text{h}}$
depends on the choice of $F$. Note that this distribution determines the
probability of finding $t_{0}$ on an ensemble of runs, in contrast with the
distributions $f$ and $f_{\theta}$ that determine the probability of opening a
brine channel within a single run. Because of this intrinsic variability, even
runs with identical bulk parameters can end up with different pond coverage.
In our simulations this variability accounted for about $5\%$ of the
variability in the minimum pond coverage, $p_{\text{min}}$.

Finally, we show that it is mainly the tail of the hole opening distribution,
$F$, that controls pond evolution. We do not know $F$ for real ice, and it can
in principle have any shape as long as it is a function monotonically
increasing from 0 to 1. The universal function, $g(\eta)$, behaves as
$g(\eta)\rightarrow 1$ for $\eta\rightarrow 0$ and $g(\eta)\rightarrow 0$ for
$\eta\rightarrow\infty$. In particular, $g$ falls below $0.1$ for $\eta\sim
10$. Therefore, when $\eta$ exceeds a value of around 10, $g$ falls to
approximately zero, and further increase in $\eta$ does not affect pond
evolution. Recall that $\eta=\eta_{0}F(\frac{t-t_{0}}{T_{\text{h}}})$, where
$\eta_{0}\equiv c\frac{l_{0}^{2}}{L^{2}}N_{0}$. If, we assume reasonable
parameter values $c\sim 3$, $l_{0}\sim 5\text{ m}$ [15], and
$\frac{N_{0}}{L^{2}}\sim 100\text{ m}^{-2}$ [6], we get an order of magnitude
estimate $\eta_{0}\sim 10000$. Therefore, $\eta\sim 10$ when $F\sim 0.001$ and
pond evolution only occurs in the tail of the distribution $F$. Thus, due to
the large density of brine channels, the exact shape of $F$ does not matter
apart from determining the weight of the distribution that falls in the tail.
We discuss this in more detail below.

Let us look at some particular examples of $F$ to see exactly how it affects
pond evolution. As a first example, let us assume that the probability density
function, $f$, falls off exponentially, $f(x)\propto e^{x}$ for $x\ll-1$. In
this case, $F(x)=Ae^{x}$ for $x\ll-1$ and $F^{-1}(y)=\ln{y/A}$ for $0<y\ll 1$,
where $A$ is a normalizing constant that also takes into account the rest of
the distribution that does not fall in the tail. If we use the approximate
relationship $t_{0}\approx-T_{\text{h}}F^{-1}(1/N_{0})$ (see Eq. S22), we find

$\eta(t)\approx
c\frac{l_{0}^{2}}{L^{2}}N_{0}Ae^{t/T_{\text{h}}-\ln{AN_{0}}}=c\frac{l_{0}^{2}}{L^{2}}e^{t/T_{\text{h}}}\text{
, }$ (S24)

for $t\ll T_{\text{h}}\ln{AN_{0}}$. Therefore, if the tail of the hole opening
distribution is exponential, pond evolution is independent of the brine
channel density, $N_{0}/L^{2}$. Pond evolution is also independent of the
shape of the distribution beyond the tail as seen from the fact that the
normalizing constant, $A$, does not appear in the equation for $\eta$. If we
recall that pond coverage during stage II is $p=p_{c}g(\eta(t))$ and note that
$g$ is approximately logarithmic in $\eta$ within the region where pond
coverage varies from $p_{c}$ to 0, we see that pond coverage falls
approximately linearly with time.

As a second example, we consider $f$ consistent with a normal distribution,
$f(x)\propto e^{-x^{2}}$, which we used in the main text to compare to
observations. In this case, $F(x)=Ae^{-x^{2}}$ for $x\ll-1$ and
$F^{-1}(y)=-\sqrt{\ln{A/y}}$ for $0<y\ll 1$, with $A$ being a normalizing
constant that accounts for the distribution beyond the tail. Following similar
logic as above, we can derive

$\eta(t)\approx
c\frac{l_{0}^{2}}{L^{2}}e^{-(t/T_{\text{h}})^{2}+2t\sqrt{\ln{AN_{0}}}/T_{\text{h}}}\text{
. }$ (S25)

Here, we see that due to the term $\sqrt{\ln{AN_{0}}}$, dependence on $N_{0}$
and $A$ is not lost. We note two facts about this dependence - 1) pond
evolution depends not on the density of brine channels but explicitly on the
number of brine channels within a drainage basin, and 2) the shape of the
distribution beyond the tail only acts to modify the effective number of brine
channels through $AN_{0}$. This unexpected dependence on the total number of
brine channels is a consequence of the fact that stage II begins when the
first hole opens rather than when a certain fraction of holes open. We
discussed this in section S8. We add that the dependence on $N_{0}$ and $A$ is
relatively weak - a 10-fold change in the effective number of brine channels
leads to only around a $10\%$ change in $\sqrt{\ln{AN_{0}}}$, and a
correspondingly small change in the rate of change of pond coverage. As noted
in Golden [6], typical brine channel densities are between $60\text{ m}^{-2}$
and $120\text{ m}^{-2}$, so the uncertainty in the number of channels is
likely not a significant source of error so long as we can reliably estimate
the size of the drainage basin, $L$. We also believe that the discussion above
partially explains why our rather ad-hoc choice of $F$ as a cumulative normal
distribution gave good agreement with observations. Namely, even if the
majority of this distribution is not normal, the prediction about pond
evolution will be accurate so long as the actual hole opening distribution
falls off as $e^{-x^{2}}$.

For the more general case of $F=Ae^{-|x|^{a}}$ in the tail, we can follow
similar logic to find that as the power $a$ increases, the dependence on
$N_{0}$ and $A$ becomes stronger. Pond coverage depends most strongly on
$N_{0}$ when the distribution $F$ has a sharp cutoff.

## References

  * Bitz and Lipscomb [1999] Cecilia M Bitz and William H Lipscomb.  An energy-conserving thermodynamic model of sea ice.  _Journal of Geophysical Research: Oceans_ , 104(C7):15669–15677, 1999.
  * Bowen et al. [2018] Brady Bowen, Courtenay Strong, and Kenneth M Golden.  Modeling the fractal geometry of arctic melt ponds using the level sets of random surfaces.  _Journal of Fractal Geometry_ , 5(2):121–142, 2018.
  * Coles [2001] Stuart Coles.  _An introduction to statistical modeling of extreme values_ , volume 208.  Springer, 2001.
  * Eicken et al. [2002] H Eicken, HR Krouse, D Kadko, and DK Perovich.  Tracer studies of pathways and rates of meltwater transport through arctic summer sea ice.  _Journal of Geophysical Research: Oceans_ , 107(C10):SHE–22, 2002.
  * Eicken et al. [2004] H Eicken, TC Grenfell, DK Perovich, JA Richter-Menge, and K Frey.  Hydraulic controls of summer arctic pack ice albedo.  _Journal of Geophysical Research: Oceans_ , 109(C8), 2004\.
  * Golden [2001] KM Golden.  Brine percolation and the transport properties of sea ice.  _Annals of glaciology_ , 33:28–36, 2001.
  * Isichenko [1992] Michael B Isichenko.  Percolation, statistical topography, and transport in random media.  _Reviews of modern physics_ , 64(4):961, 1992\.
  * Landy et al. [2014] Jack Landy, Jens Ehn, Megan Shields, and David Barber.  Surface and melt pond evolution on landfast first-year sea ice in the canadian arctic archipelago.  _Journal of Geophysical Research: Oceans_ , 119(5):3054–3075, 2014.
  * Lüthje et al. [2006] Mikael Lüthje, DL Feltham, PD Taylor, and MG Worster.  Modeling the summertime evolution of sea-ice melt ponds.  _Journal of Geophysical Research: Oceans_ , 111(C2), 2006\.  ISSN 2156-2202.  doi: 10.1029/2004JC002818.  C02001.
  * Ono [1967] Nobuo Ono.  Specific heat and heat of fusion of sea ice.  _Physics of Snow and Ice_ , 1(1):599–610, 1967\.
  * Polashenski et al. [2012] Chris Polashenski, Donald Perovich, and Zoe Courville.  The mechanisms of sea ice melt pond formation and evolution.  _Journal of Geophysical Research: Oceans_ , 117(C1), 2012\.
  * Polashenski et al. [2017] Chris Polashenski, Kenneth M Golden, Donald K Perovich, Eric Skyllingstad, Alexandra Arnsten, Carolyn Stwertka, and Nicholas Wright.  Percolation blockage: A process that enables melt pond formation on first year arctic sea ice.  _Journal of Geophysical Research: Oceans_ , 122(1):413–440, 2017.
  * Popović and Abbot [2017] Predrag Popović and Dorian Abbot.  A simple model for the evolution of melt pond coverage on permeable arctic sea ice.  _The Cryosphere_ , 11(3):1149–1172, 2017.
  * Popović et al. [2018] Predrag Popović, BB Cael, Mary Silber, and Dorian S Abbot.  Simple rules govern the patterns of arctic sea ice melt ponds.  _Physical review letters_ , 120(14):148701, 2018\.
  * Popović et al. [2020] Predrag Popović, Justin Finkel, Mary C Silber, and Dorian S Abbot.  Snow topography on undeformed arctic sea ice captured by an idealized “snow dune” model.  _Journal of Geophysical Research: Oceans_ , 125(9):e2019JC016034, 2020.
  * Rösel et al. [2015] Anja Rösel, Lars Kaleschke, and Stefan Kern.  Gridded melt pond cover fraction on arctic sea ice derived from terra-modis 8-day composite reflectance data bias corrected version 02, 2015.
  * Scagliarini et al. [2018] Andrea Scagliarini, Enrico Calzavarini, Daniela Mansutti, and Federico Toschi.  Modelling sea ice and melt ponds evolution: sensitivity to microscale heat transfer mechanisms.  _arXiv preprint arXiv:1809.06924_ , 2018.
  * Skyllingstad et al. [2015] Eric D Skyllingstad, Karen M Shell, Lee Collins, and Chris Polashenski.  Simulation of the melt season using a resolved sea ice model with snow cover and melt ponds.  _Journal of Geophysical Research: Oceans_ , 120(7):5194–5215, 2015.
  * Untersteiner [1961] N Untersteiner.  On the mass and heat budget of arctic sea ice.  _Archiv für Meteorologie, Geophysik und Bioklimatologie, Serie A_ , 12(2):151–182, 1961.
  * Wakatsuchi and Saito [1985] Masaaki Wakatsuchi and Takashi Saito.  On brine drainage channels of young sea ice.  _Annals of glaciology_ , 6:200–202, 1985.