A simple model for the evolution of melt pond coverage on permeable Arctic sea ice

As the melt season progresses, sea ice in the Arctic often becomes permeable enough to allow for nearly complete drainage of meltwater that has collected on the ice surface. Melt ponds that remain after drainage are hydraulically connected to the ocean and correspond to regions of sea ice whose surface is below sea level. We present a simple model for the evolution of melt pond coverage on 5 such permeable sea ice floes in which we allow for spatially varying ice melt rates and assume the whole floe is in hydrostatic balance. The model is represented by two simple ordinary differential equations, where the rate of change of pond coverage depends on the pond coverage. All the physical parameters of the system are summarized by four strengths that control the relative importance of the terms in the equations. The model both fits observations and allows us to understand the behavior 10 of melt ponds in a way that is often not possible with more complex models. Examples of insights we can gain from the model are: (1) the pond growth rate is more sensitive to changes in bare sea ice albedo than changes in pond albedo, (2) ponds grow slower on smoother ice, and (3) ponds respond strongest to freeboard sinking on first year ice and sidewall melting on multiyear ice. We also show that under a global warming scenario, pond coverage would increase, decreasing the overall ice 15 albedo, and leading to ice thinning that is likely comparable to thinning due to direct forcing. Since melt pond coverage is one of the key parameters controlling the albedo of sea ice, understanding the mechanisms that control the distribution of pond coverage will help improve large-scale model parameterizations and sea ice forecasts in a warming climate.


Introduction
Over the past 40 years, Arctic summer sea ice extent has reduced by 50 %, making it one of the most sensitive indicators of man-made climate change (Serreze and Stroeve, 2015;Stroeve et al., 2007;Perovich and Richter-Menge, 2009).This rapid decrease is at least partially due to the ice-albedo feedback (Zhang et al., 2008;Screen and Simmonds, 2010;Perovich et al., 2007).Moreover, if the ice-albedo feedback is strong enough it could lead to instabilities and abrupt changes in ice coverage in the future (North, 1984;Holland et al., 2006;Eisenman and Wettlaufer, 2008;Abbot et al., 2011).The albedo of ice is significantly reduced by the presence of melt ponds on its surface (Eicken et al., 2004;Perovich and Polashenski, 2012;Yackel et al., 2000).Therefore, understanding the evolution of melt ponds is essential for understanding the ice-albedo feedback and, consequently, the evolution of Arctic sea ice cover in a warming world.This means that accurate melt pond parameterizations must be incorporated into global climate models (GCMs) to improve their sea ice forecasts (Flocco et al., 2010;Holland et al., 2012;Pedersen et al., 2009).The main difficulties with including accurate melt pond parameterizations in large-scale models are that pond evolution is nonlinear and that it is the result of a variety of different physical processes operating on a range of length and time scales.For these reasons, it is important to understand the mechanisms that drive the evolution of melt ponds.
Typically, the evolution of pond coverage on first-year ice proceeds in fairly consistent stages (Polashenski et al., 2012;Perovich et al., 2003;Landy et al., 2014;Webster et al., 2015).First the ponds grow quickly while the ice is impermeable.Next they drain quickly and pond coverage shrinks as the ice transitions from impermeable to permeable.Then the ponds grow slowly while the ice is permeable and pond water remains at sea level.Finally, the ponds either refreeze or the floe breaks up.The stage when ice is highly permeable is typically the longest, often longer than the first two stages combined.This stage is particularly suitable to model, since the ponds can be assumed to be at sea level and hydraulically connected to the ocean.On multiyear ice, ponds also experience a growth and a drainage stage, but often do not drain to sea level.On some occasions, however, ponds on multiyear ice can drain to sea level as well.
In this paper we will present a simple "0-D" model for the evolution of melt pond coverage on sea ice floes.We will assume that ice is permeable, ponds are at sea level and hydraulically connected to the ocean, the whole ice floe is in hydrostatic balance, and different points on the ice surface may melt at different rates.The purpose of our model is (1) to clarify the roles in the evolution of pond coverage played by energy fluxes, the ice thickness, bulk ice density, ice roughness, and initial pond coverage; (2) to provide a simple, yet accurate, way to estimate the pond coverage as a function of time; (3) to understand the behavior of melt ponds under general environmental conditions; and (4) to investigate different types of qualitative behavior that can arise from differential melting and maintaining hydrostatic balance.Skyllingstad et al. (2009) also describe pond growth on permeable ice, but they include only pond growth by lateral melt of pond walls.This contrasts with our model, which includes pond growth by vertical changes of the topography.Our models are different, but complementary, and we will draw parallels between our two models when discussing the possibility of lateral melt.Aside from Skyllingstad et al. (2009), previous melt pond modeling efforts include works by Taylor and Feltham (2004), Lüthje et al. (2006), Scott and Feltham (2010), and Flocco and Feltham (2007), who all created comprehensive models that allowed for more realistic representations of physical processes such as heat and salt balance, and meltwater routing and drainage.The advantage of our model is its simplicity, which makes it possible to clarify the roles of each of the physical parameters involved.
This paper is organized in the following way.In Sect. 2 we build a simple model for the evolution of pond coverage.In Sect.3, we compare the model to observations.In Sect. 4 we discuss realistic values of physical parameters and solve the model numerically.In Sect. 5 we assess the impacts of sea ice roughness and develop a simple parameterization to estimate mean pond coverage after a certain amount of time without solving the model.In Sect.6 we analyze the model analytically to gain a better understanding of the factors influencing pond evolution.In Sect.7 we discuss lateral melt and internal melt combined with effect of density variations.Finally in Sect.8 we summarize our results and conclude.In Appendices A, B, C, and D we discuss some of the more technical aspects of our model.

Building the simple 0-D model
In this section, we build the model for the evolution of melt pond coverage and then solve it using realistic physical parameters.Before we proceed to build the quantitative model, we will first state the assumptions and discuss the physical mechanisms driving pond evolution.

Assumptions of the model
Our model focuses on the stage of pond evolution when ice is highly permeable and all the meltwater created can be quickly removed to the ocean.The beginning of this stage can be identified as the point in time when the meltwater on the ice surface has drained to sea level, such that the remaining ponds correspond to places on the ice surface that are below sea level.We will assume that from this point on, the ponds are hydraulically connected with the ocean, and the only way for pond coverage to increase is for the points on the ice surface which were above sea level to sink or melt below sea level.In reality, ponds can also grow through horizontal melting of their sidewalls.As some observations suggest that this type of growth is small at least on first-year ice (Polashenski et al., 2012;Landy et al., 2014), we neglect it (see Sect. 7.1 for further discussion).Furthermore, we will assume that all the melt occurs at the surface or the bottom of the ice.We thereby neglect the possibility of internal melt.We will also assume that ice has a uniform bulk density throughout the vertical column, and we discuss the effects of vertical nonuniformity in bulk density together with effects of internal melt in Sect.7.2.Finally, we will assume that the entire ice floe is in hydrostatic balance.
The main goal of our model is to determine the fraction of the ice surface above sea level that falls below sea level after some time.Therefore, we focus on the vertical displacements of points on the surface of the ice in response to melt.To this end, we define the ice topography, s(r), as the elevation of the ice surface above sea level at the point r, and we define melt ponds as those regions where s(r) < 0. There are two main reasons why the topography might change in response to ice melt: 1. First, the topography at a point r at the surface changes when ice at that point melts (Fig. 1a).Here, the rate of change of topography at a point depends only on local characteristics of that particular point.For this reason, we will call this type of motion "local."Points on the surface that melt locally move "downwards," i.e., to lower elevations above sea level.
2. Second, in order to maintain hydrostatic balance, the entire ice surface can shift up or down in response to mass being removed above or below sea level.Since we are assuming that the entire ice floe is in hydrostatic balance, melting any region of ice moves the entire floe as a rigid body (Fig. 1b).For this reason, we will call this type of motion the "rigid body" motion.Melting above sea level induces an upward rigid body motion, whereas melting below sea level induces a downward rigid body motion.An ice floe is not a rigid body, but up to its flexural wavelength (roughly 30 m on 1.5 m thick ice) we can approximate it as such.As the flexural wavelength is larger than the typical scale of melt ponds (roughly 10 m), the rigid body approximation is likely good.
At each point on the ice surface the change in elevation above sea level can be calculated as the sum of these two contributions.
In our model, ponds grow in two ways, "freeboard sinking" and "enhanced melting": 1. Freeboard sinking represents the average change in freeboard height (average height above sea level of bare ice).In this way the topography of ice above sea level remains unchanged.Freeboard sinking should not be confused with rigid body motion: the average freeboard height always decreases as a response to ice thinning, whereas the rigid body motion can point both upward and downward depending on whether mass is lost above or below sea level.Both rigid body motion and average local melting contribute to freeboard sinking.
2. Enhanced melting represents the change in the shape of the topography without changing its average height.Ponds can grow in this way if some regions melt faster than average.Therefore, a positive deviation in the local melt rate can grow ponds.Conversely, a negative deviation in the local melt rate can slow down or even reverse pond growth.Pond growth only occurs due to topography changes near sea level.Therefore, deviations from the mean melt rate for points high above the sea level do not influence pond evolution since these points are correlated with points close to sea level only through hydrostatic adjustment, which is determined by the average melt rates rather than the deviations from the average.

Equation for the evolution of topography
We now proceed to build the quantitative model of pond evolution.Following the above ideas, we divide the total rate of change of vertical position of the point r on the surface of the ice, ds dt (r), into a contribution from rigid body motion, ds rigid body dt , and a contribution from local melting, ds loc dt (r): Ice above sea level (asl) must hydrostatically balance ice below sea level (bsl).We can write this hydrostatic balance as where m asl and m bsl represent the mass of ice above and below sea level, and ρ w and ρ i represent the densities of sea water and pure ice.Throughout the paper we use ρ w = 1025 kg m −3 and ρ i = 916 kg m −3 .The mass above and below sea level can change either because the ice melts or because the floe moves as a rigid body, changing the proportion of ice above and below sea level.Therefore, differentiating Eq. ( 2) and splitting into melt and rigid body contributions, we find dm where l = 334 kJ kg −1 is the latent heat of melting, F bi is the total energy flux used for melting bare ice averaged over all bare ice, F mp is the total energy flux used for melting ponded ice averaged over ponded ice, and F bot is the total energy flux used for melting the ice bottom averaged over the ice bottom.A bi , A mp , and A are the area of bare ice, the area of melt ponds, and the area of the entire floe, respectively.Since floating up or down does not change the total mass of the ice, mass changes above and below sea level due to rigid body motion are equal with an opposite sign, dm where ρ b is the bulk ice density.This is the density of sea ice once all the brine has drained and is always less than ρ i .
We assume it to be uniform throughout the vertical ice column, but we discuss the effects of vertical variations in ρ b in Sect.7.2.Substituting Eqs. ( 4) and ( 5) into Eq.( 3), solving for ds rigid body , and differentiating with respect to time, we find the rate of change of surface topography due to rigid body motion to be ds rigid body dt = The three terms in large square brackets correspond to topography change due to bare ice melting, ponded ice melting, and ice bottom melting.Rigid body motion depends only on spatially averaged energy fluxes, which in turn depend on parameters such as the average insolation on the floe, the average albedo, and the average longwave, sensible, latent, and bottom heat fluxes.If bare and ponded ice melt only from energy absorbed by the upper surface of the ice, the fluxes F bi and F mp can also be written in terms of albedo as where α bi and α mp are the average albedos of bare and ponded ice, F sol is the solar flux, and F r is equal to the sum of net longwave, net sensible, and net latent heat fluxes.This parameterization neglects light transmission and assumes that all of the energy is deposited in the surface.Much of the variation in albedo of ponded ice is due to the fact that the pond bottom is partially transparent, and energy is deposited in the ocean instead of directly in the ice.However, this does not make much difference in our model since the energy deposited in the ocean is likely used for melting ice below sea level anyway.
Local displacement, ds loc , quantifies how much the ice surface topography changes as a result of local melt.We can determine the local melt rate from F surf (r), the flux of energy used for melting the ice surface at a point r: where the positive direction is defined as upwards.The local flux depends on parameters such as the local albedo, the local insolation, the local longwave, sensible and latent heat fluxes, and the angle between ice and incoming radiation at that point.
The flux F surf (r) averaged over all the points on the surface of the ice above sea level equals F bi : where < . . .> represents averaging over all the points on bare ice.For this reason, we will parameterize the rate of local melting as where k(r) is a nondimensional number that quantifies the deviation of the melt rate at the point r from the mean melt rate of the bare ice surface, which depends on the detailed conditions of ice and its environment.The parameter k could be either greater than or less than 1.Here we will take k to be constant in time, but in reality it need not be.Finally, according to Eq. (1) we add Eqs. ( 6) and (11) to get the equation for the evolution of the bare ice topography.We express this in terms of melt pond fraction, x ≡ Here, we split the equation into two terms, enclosed by the square brackets.The first term represents the local deviation from the average surface melt rate, which changes the general shape of the topography while preserving its average height above sea level.We identify this term with enhanced melting.The second term represents a global shift of the average elevation above sea level due to freeboard sinking.
In this way, the topographic evolution equation can be split into two terms, enhanced melting and freeboard sinking: where ds em dt and ds fs dt are contributions from enhanced melting and freeboard sinking, and they correspond to the first and second term of Eq. (12).

Model for the evolution of pond coverage
We now need to relate the vertical displacements near the sea level to the change in area of the melt ponds.To this end we define the hypsographic curve, s(x h ), which relates the elevation above sea level, s, to the percent of ice surface below that elevation, x h (Fig. 2).Such curves have been measured and reported on several occasions (e.g., Fig. 8 of Eicken et al., 2004, or Fig. 8 of Landy et al., 2014).If the ice is highly permeable, the melt pond fraction, x, can be inferred from a hypsographic curve as the intersection of sea level with the curve.Since ponds are hydraulically connected with the ocean, the average freeboard height of bare ice, h, depends on the pond fraction.The average freeboard height, h, can be expressed in terms of the ice thickness H and the pond fraction as Here, the average freeboard height is defined as the elevation of the ice surface above sea level averaged over bare ice.For two ice floes of the same thickness, the one with higher pond coverage will also need to have a higher average freeboard in order to maintain hydrostatic balance.
The above sea level part of every measured hypsographic curve we tested can be fit relatively well with a tangent function (Fig. 2a, red line).We will assume that this fit holds for a wide range of different sea ice floes and use it to initialize our model with different physical parameters.We give the exact form of this function in Appendix A (Eq. A1).To get a hypsographic curve for a particular initial pond fraction, x i , and ice thickness, H , we set it to zero at the initial pond coverage, s(x h = x i ) = 0, and rescale it vertically to get a freeboard that hydrostatically balances the floe.The topography below sea level is not important for the evolution of pond coverage if the pond coverage grows, and we replace it with a straight line.
We show several curves for different initial ice thickness and initial pond coverage in Fig. 2b and c.We note that the initial pond fraction, x i , corresponds to the pond fraction when ice first becomes permeable.Once we choose x i and H , the tangent function Eq. (A1) has only two unconstrained parameters, p 1 and p 2 , that determine the exact shape of the curve.Knowing additional physical parameters, such as ice roughness, we can constrain additional parameters of this curve.Throughout this paper we will mostly use p 1 and p 2 that fit the measurements of the hypsographic curve made by Landy et al. (2014) for 25 June 2011 or the measurements made during the SHEBA mission along the topography profile "1" on 10 July 1998.However, when examining the effects of sea ice roughness, we will vary these parameters to get curves of different shape.Several examples of hypsographic with different p 1 and p 2 are shown in Fig. 2d.
In the case of pure freeboard sinking the overall shape of the hypsographic curve does not change as the ice melts.In- stead the whole curve is shifted following a displacement of ds fs (Fig. 3a).We can calculate the resulting change in pond coverage as where ds fs is the vertical displacement of the bare ice topography due to freeboard sinking (as determined by the second term in Eq. 12), and dx h ds (x) is the change in pond fraction for a vertical shift of the ice surface of ds fs when the pond fraction is equal to x.It is equal to the reciprocal of the derivative of the hypsographic curve, s(x h ), evaluated at x h = x.Substituting ds fs dt from Eq. ( 12) we find dx dt = where x ≡ x x i and 1 − x ≡ 1−x 1−x i are the pond and bare ice fractions normalized by the initial pond and bare ice fractions, and  of the hypsographic curve.We have defined the strengths of pond growth by freeboard sinking due to melting bare, ponded, and ice bottom, S bi , S mp , and S bot , as The nondimensional factors x, 1 − x, and d xh dŝ (x) are chosen to be of the order unity, so that S bi , S mp , and S bot control the strengths of pond growth by melting bare ice, melting ponded ice, and melting ice bottom.The reciprocals of the strengths represent the timescales of the growth modes.
The set of parameters needed to describe pure freeboard sinking can be further reduced by rewriting Eq. ( 16) as where S 1 ≡ S mp −x i S bi /(1−x i ) and S 2 ≡ S bot +S bi /(1−x i ) represent a minimal set of parameters needed to describe pure freeboard sinking.However, these parameters do not have a clear physical interpretation, and we will henceforth focus only on S bi , S mp , and S bot .
Next we need to consider the contribution from enhanced melting.Before doing so we need to make some assumptions about the nature of enhanced melt.There are multiple physical processes that can cause the melt rate to deviate from the mean.One process that stands out as being particularly important is albedo decrease due to ice wetting: ice close to sea level will likely be wet and therefore have a lower albedo compared to ice higher up.The deviation from the mean melt rate in this case depends primarily on the height above sea level.Another potential contribution to height-dependent enhanced melt may effectively come from random fluctuations in the melt rate around the average: ice near the sea level has a higher probability of falling below sea level due to random fluctuations than ice higher up.After falling below sea level, ice becomes ponded, melts faster, and is unable to return to its previous position.Other processes, such as lateral melt, may not depend on height above sea level, but for now we neglect this possibility (see Sect. 7.1 for discussion).
Because of the processes described above, we will assume that the deviation from the mean melt rate, k(r) − 1, depends only on height above sea level, s.In this scenario, we need to consider enhanced melting together with freeboard sinking, as freeboard sinking constantly supplies new ice to low elevations to be affected by enhanced melting.Effects of enhanced melting and freeboard sinking can be approximately separated if, instead of height dependence, enhanced melting is constrained to act on a fixed fraction of bare ice.In this case, a constant fraction of bare ice that would experience enhanced melting would evolve, at least approximately, independently of freeboard sinking.
Therefore, we will consider two cases of enhanced melting.Firstly, we will consider a height-dependent enhanced melting.In particular, we will assume that k(0 < s < s) ≡ k and k(s > s) ≡ 1, where s is a height above which there is no enhanced melting and below which enhanced melting is constant k > 1.This is the case we ultimately wish to describe.We describe a potential model for pond growth under this assumption in Appendix B and Fig. 3d.However, from a practical viewpoint, it is simpler to consider enhanced melting which acts upon a fixed fraction of bare ice.In this case, we will assume that k(x < x h < x + δ) ≡ k and k(x h > x + δ) ≡ 1, where δ is a fraction of ice affected by enhanced melting (Fig. 3b).In Appendix B, we show that, if δ is appropriately chosen, a height-dependent model and a fixed fraction model become equivalent.Therefore, we will first solve a model assuming a fixed δ and no freeboard sinking and then relate it to a fixed s model by choosing the appropriate δ.
We note that the assumption that k(r) = 1 high above the sea level and k(r) > 1 near the sea level is strictly not true since averaged over all of bare ice k(r) needs to equal 1.However, it is approximately true if s or δ are small, such that the area where k(r) = 1 is small compared to the total area of bare ice.Also, we have assumed k(r) = 1 high above the sea level without loss of generality, since deviations from the mean melt rate high above the sea level are not important, as only ice close to sea level may become ponded.Now we proceed to consider the case of "pure enhanced melting" that assumes a fixed fraction of the ice, δ, melts, and there is no freeboard sinking (Fig. 3b).If there is no topographic variation above sea level, and the entire ice floe above sea level has the same height, h, the pond coverage would grow by δ after a time t = h ds em /dt , where ds em /dt is the rate of change of topography due to enhanced melting as determined by the first term of Eq. ( 12).Therefore, the pond growth rate in this case would be x t = δ h ds em dt .If there is non-negligible topography above sea level described by the hypsographic curve, the time t it takes for pond coverage to grow by δ would be t = s(x h =x+δ) ds em /dt .Here, s(x h = x + δ) is the original hypsographic curve evaluated at x h = x + δ.We will assume this expression generally holds for enhanced melting.Thus, we arrive at the expression for pond growth due to pure enhanced melting with fixed δ: If δ is small compared to the variation in the hypsographic curve, we can substitute s(x + δ) with s(x).This is only not justified near the beginning of the melt, when s(x) ≈ 0. Substituting ds em dt from Eq. ( 12) we find where ŝ(x) ≡ s(x) h is the nondimensional hypsographic curve, and the strength of the enhanced melting, S em , is defined as Ultimately, however, our goal was to describe the heightdependent enhanced melting.In Appendix B, we showed that such a model can be approximated with a fixed fraction model, if we appropriately relate δ and s.Here we simply state the result: Here, ds em ds fs represents the ratio of the topographic rate of change due to enhanced melting to freeboard sinking and is given by where |F | are the representative values of energy fluxes, e.g., their time averages.Therefore, the strength of heightdependent enhanced melting becomes We have made a number of assumptions in deriving the expression for enhanced melting.Below we compare this model to a more complicated "1-D" model and show that all these assumptions are justified.We also show that if the function describing the local melt rate, k(s), has a nontrivial dependence on height above sea level, parameter S em is better replaced with a parameter: In this way, we have separated the effects of freeboard sinking and enhanced melting.Finally, we will assume that contributions from freeboard sinking and enhanced melting can be added independently.Therefore, we solve Eq. ( 16) for pure freeboard sinking and Eq. ( 20) for enhanced melting independently, and we add them together to get the full evolution of pond coverage, x(t): where x fs (t) and x em (t) are solutions to Eqs. ( 16) and ( 20), both forced using the same parameters and initialized with the same initial pond fraction x i .This concludes the 0-D model.Equation ( 26) represents a sum of solutions to two simple ordinary differential equations, in which the rate of change of pond fraction depends on the pond fraction.Here, we have reduced the number of parameters from the original 10 (H , x i , ρ b , F bot , F sol , F r , α bi , α mp , k, and s) to 4 (S bi , S mp , S bot , and S em ).The strengths of freeboard sinking, S bi , S mp , and S bot , depend only on the parameters that are available in GCM simulations and are relatively easily measured in observational studies.The enhanced melting strength, S em , however, also depends on the difficult-to-measure parameters k and s that describe the melt rate near the sea level and may also have contributions from processes that are not height dependent.Furthermore, as we discuss below, ice roughness can also play an important role in pond evolution.With reliable constraints on these parameters, our model would be a useful parameterization in GCMs for pond growth after ice becomes permeable.

Testing the model
In order to test the assumptions we made to simplify the model, we have developed a "1-D" model in which we explicitly determine pond evolution when both freeboard sinking and enhanced melting are happening simultaneously.Apart from resolving the melt rates in one dimension, the underlying assumptions for the 1-D model are essentially the same as for the simple model.For this reason, we simply give an outline for this model, without discussing it in much detail.
In the 1-D model, we evolve the hypsographic curve by prescribing a melt rate, ds loc , to each point on the hypsographic curve depending on the height above sea level (Fig. 3c).The hypsographic curve high above sea level melts at a uniform rate, whereas the hypsographic curve slightly above sea level melts at an enhanced rate.Parts of the curve below sea level melt at a uniform rate determined by the flux used for melting ponded ice, F mp .Finally, hydrostatic adjustment is calculated by finding the ice thickness directly at each time step and placing the floe in hydrostatic balance.The evolution of pond coverage obtained from this model is shown in Fig. 4a.The comparison with the simple 0-D model is excellent.
The 1-D model allows us some freedom to test the detailed assumptions of the 0-D model.First, we can test how the functional form of k(s) affects the pond evolution (Fig. 4b).
The functions k(s) were chosen such that they all have the same integral parameter < S em > defined in Eq. ( 25). Figure 4a shows that in each of these cases the evolution of pond coverage proceeds nearly identically.Second, we can test the difference between an assumption that enhanced melting acts below a constant height s and an assumption that enhanced melting acts on a constant fraction of ice, δ.The yellow line in Fig. 4a shows that if δ and s are chosen according to Eq. ( 22), both assumptions yield very similar results.Finally, we can test the effects of varying pond albedo.In reality pond albedo decreases as the ponds deepen.We assume a dependence of pond albedo on pond depth reported in Table VII of Morassutti and Ledrew (1996) for mean broadband albedo.The magenta line in Fig. 4a shows that allowing for pond albedo to vary has a negligible effect on pond evolution.
We should note that, when both freeboard sinking and enhanced melting occur simultaneously, the agreement between the 0-D model and the 1-D model becomes poor if the hypsographic curve is convex (e.g., Fig. 2d, blue curve), and the 0-D model should be used with care.Happily, the measured hypsographic curves are mostly concave, in which case the agreement between the two models is excellent.In Fig. 5, we compare the results from our model to observations made on a 200 m long albedo line during SHEBA (red line).Ice along the albedo line was level multiyear ice, but the ponds drained to sea level after some time, which makes them amenable to our model (Perovich et al., 2003).
The pond coverage along the albedo line dropped to a minimum around the end of June.Therefore, we choose to model only the period after 1 July.In order to keep the albedo line pristine, no thickness measurements were made.However, relatively close to the albedo line, topography measurements were made along a level multiyear ice profile roughly every 10 days.After approximately 10 July, ponds along the topography profile also drained to sea level.We show the topography profile pond coverage in blue dots (we have artificially subtracted 0.05 from the pond coverage to facilitate comparison with the pond coverage along the albedo line).
The pond coverage both along the topography profile and along the albedo line follows roughly the same trend, suggesting that the physical parameters driving the pond evolution in the two places are likely similar.Based on the average freeboard height, we estimate the ice thickness on 10 July to be roughly 1.4 m along the topography profile, meaning that on 1 July ice thickness was around 1.6 m.Therefore, we assume the same thickness for the ice along the albedo line and use a hypsographic curve corresponding to the one measured along the topography profile on 10 July (Fig. 2a, dashed line).In order to run our model, we use the melt rates of bare ice, ponded ice, and ice bottom measured directly using ablation stakes during SHEBA (Perovich et al., 2003).
We choose a realistic ρ b = 850 kg m −3 (Timco and Frederking, 1996).We have no way of directly constraining the parameters s and k that control the strength of enhanced melting.Therefore, we treat S em as a fitting parameter.Choosing S em = 0.22 month −1 fits the observations well by eye.This value can be obtained using s = 15 cm and k = 1.7, which likely fall at the upper end of the range of reasonable values for these constants (see Sect. 4 for a discussion on s and k).Such a high value of S em can be explained by a significant contribution from lateral melting.The full black line in Fig. 5 represents a solution to the full Eq.( 26).The agreement between model and observation is excellent, with a maximum discrepancy of 3 % pond coverage at the end of the melt season.The dashed black line represents the contribution to pond growth due to freeboard sinking, whereas the dotted line corresponds to enhanced melting.Almost all pond growth in this case is due to enhanced melting.This is due to ice topography.On multiyear ice, meltwater typically collects in depressions formed by ponds in previous years.The topography created in this way is highly bimodal, and, after drainage, ponds typically have steep walls.Bare ice topography, in contrast, is relatively smooth, preventing new pond formation.This is ap-parent in the hypsographic curve we used.Such a topography inhibits freeboard sinking, and pond coverage grows mostly by enhanced melting acting near the pond sidewalls, growing the existing ponds.In addition to height-dependent enhanced melting we introduced in the previous section, in this case there is likely a significant contribution from lateral melting as well.This contribution helps explain the high value of S em we had to choose to get a close agreement between our model and observations.First-year ice topography, in contrast, permits ample pond growth through freeboard sinking.Observations suggest that on first-year ice ponds grow primarily due to freeboard sinking (Polashenski et al., 2012;Landy et al., 2014).

Numerical solutions
We now solve Eq. ( 26) numerically to gain intuition about the behavior of our model.We use a set of realistic parameters we will henceforth refer to as the "default parameters." For shortwave, longwave, latent, and sensible heat fluxes, we use values inferred by Skyllingstad et al. (2009) using hourly measurements from the SHEBA mission.We use the bottom heat flux inferred from measurements of ice bottom ablation during the SHEBA mission (Perovich et al., 2003).The albedo of bare ice can vary between 0.5 and 0.7 (Hanesiak et al., 2001), while the albedo of melt ponds can vary between 0.1 and 0.6, depending on pond depth and conditions of ice at the pond bottom (Morassutti and Ledrew, 1996;Perovich et al., 1998;Perovich, 1996).Here we prescribe a default bare ice albedo of 0.55 and a default pond albedo of 0.2.We use a realistic bulk ice density of ρ b = 850 kg m −3 (Timco and Frederking, 1996).We use an initial ice thickness of 1.5 m and use the first-year ice topography measured by Landy et al. (2014) adjusted for the prescribed ice thickness and initial pond fraction (usually x i = 0.2).We will assume enhanced melting is entirely due the albedo dependence on height above sea level.Some preliminary results based on field measurements of bare ice albedo on first-year ice suggest that albedo changes from around 0.3 near sea level to around 0.55 at a height of around 10 cm above sea level, after which the correlation between albedo and surface elevation tapers off (Chris Polashenski, personal communication, 2017).Using such an albedo and the average values of shortwave, longwave, latent, and sensible heat fluxes, we can estimate the rate of melt as a function of height above sea level, k(s) = F (s) F bi .Using Eq. ( 25), we can then find the integral parameter < S em >.We choose s = 6 cm and k = 1.7 to correspond to the same integral parameter.We should note that there is significant scatter in the data, and measurements correspond to only one study.Therefore, this is a rough estimate of enhanced melting, but it is likely of the correct order of magnitude.
Figure 6a shows the solution to Eq. ( 26) for different initial conditions.We can see that ponds grow more rapidly when the initial pond coverage is lower, and the pond evolution curves cluster together as time progresses.This is because lower initial pond coverage corresponds to lower initial freeboard height, making the pond growth more rapid.The dashed line corresponds to the solution using the fluxes timeaveraged over the 30-day run.The solutions using the averaged fluxes are very similar to the ones using time-varying fluxes, meaning that daily and even monthly variations in the forcing have little effect on pond growth.This insensitivity to short timescale variations in the forcing means that pond coverage evolution may be faithfully represented in the large-scale models, as it would not be affected by the coarse timescales of those models.Henceforth, we will use the timeaveraged fluxes.
A larger ice thickness means a higher freeboard.For this reason, ponds grow more slowly on thicker ice.Because the pond growth rate is inversely proportional to ice thickness, pond coverage is more sensitive to variations in ice thickness when the ice is thin (Fig. 6b).In Fig. 6b we see that a 0.5 m difference in the initial ice thickness (between a floe 1.5 m and a floe 2 m thick) can mean a 20 % difference in pond coverage at the end of the melt season.
Figure 6c shows the dependence of pond coverage on albedo.A variation of 0.1 in bare ice albedo has a much larger effect on pond evolution than the same change in pond albedo.The reason is that melting ponded ice only affects pond coverage through downward rigid body motion of the floe, whereas melting bare ice grows the ponds through both enhanced melting and freeboard sinking.Furthermore, when pond coverage is low, rigid body motion due to ponded ice melting is less efficient than that due to bare ice melting because it is proportional to melt pond fraction.
The parameters controlling the strength of enhanced melting are the least constrained parameters in our model.In Fig. 6d we show the dependence of pond evolution on the height below which enhanced melting is active, s.Exploring a range of realistic values for s, 0 < s < 15 cm, we find that the pond fraction at the end of the melt season can vary by about 30 %.This difference would be larger if we chose a smaller ice thickness.The effects of changing k are relatively small, so long as k is large enough (not shown).For example, using current parameters, pond coverage evolution becomes fairly insensitive to k when k > 1.5.Smaller values of k, however, can significantly impact pond evolution.If k is sufficiently smaller than 1, S em can become negative and the pond coverage can stop growing.In this case, ice near the sea level melts slowly enough such that an upward rigid body movement due to melting ice high above sea level pushes the ice near sea level upwards, preventing pond coverage growth.The evolution of such a pond coverage cannot be represented well in our model since the equation for enhanced melting becomes invalid in this case, and the blue curve in Fig. 6d serves therefore simply as an illustration.

Pond evolution is slower on smoother ice
The evolution of pond coverage in our model depends on the detailed shape of the hypsographic curve, which is not captured by the strengths of freeboard sinking and enhanced melting.As we show below, pond coverage is sensitive to such details and in particular to ice roughness.Below we will introduce the "effective strengths", S * , which approximately capture the effects of roughness and allow us to estimate mean pond coverage after a period of time.Using effective strengths, we will demonstrate how multiyear ice topography suppresses pond growth by freeboard sinking, while first-year ice topography permits it.
In the tangent function parameterization, Eq. (A1), the exact shape of the hypsographic curve is determined by parameters p 1 and p 2 .Here, we will not discuss these parameters individually but will rather focus on often measured bare ice roughness, σ , defined as the standard deviation of surface elevation of ice above sea level: We will use the nondimensional form of bare ice roughness, defined as σ ≡ σ h .Typically, a concave hypsographic curve (e.g., Fig. 2d, red curve) will have a small σ , whereas a convex hypsographic curve (e.g., Fig. 2d, blue curve) will have a high σ .
During the permeable stage, all else equal, ponds will grow more rapidly on rougher ice, since a larger fraction of ice is close to sea level.This is not true on impermeable ice, as meltwater filling deep topographic lows on rough ice will cover a smaller area relative to the same amount of meltwater filling shallow topographic lows on smooth ice.For this reason, the initial pond coverage will likely be smaller on rougher ice due to a smaller pond coverage during the impermeable stage.
Figure 7 shows the pond coverage evolution due to pure freeboard sinking (Fig. 7a) and pure enhanced melting (Fig. 7b) for hypsographic curves with different parameters p 1 and p 2 and all other parameters kept constant.For each choice of p 1 and p 2 , we find the normalized bare ice roughness, σ , represented by the color of the curves.Blue colors correspond to low roughness and red colors to high roughness.Pond evolution on measured topographies (Fig. 2a) is also shown.We can see that although roughness does not fully determine the pond evolution, it is a viable proxy for how pond coverage will evolve, with high roughness curves typically having a higher average pond coverage.We wish to quantify the effect of roughness by its impact on the mean pond coverage.In particular, we hope to find the "effective strengths", S * ( σ ), which include the roughness effects and allow us to easily estimate the average pond coverage after some time t: where < x(t) >≡ t 0 x(t)dt t . Effective strengths are proportional to strengths of freeboard sinking and enhanced melting we derived in Sect.2.3.In general they themselves may depend on time and are independent of time only if pond coverage evolution is linear, x(t) = St + x i , in which case S * = S, where S is either S fs ≡ (S bi + S mp + S bot ) in the case of freeboard sinking or S em in case of enhanced melting.
In Appendix C, we describe the procedure to estimate the effective strengths as functions of nondimensional roughness and time.Here, we only state the result: where S * fs is the effective strength of freeboard sinking, S * em is the effective strength of enhanced melting, and tem ≡ S em t 1−x i is the nondimensional time of pond evolution due to enhanced melting.The terms in square brackets represent the corrections due to roughness.If both freeboard sinking and enhanced melting occur simultaneously the total effective strength is the sum of these two, S * = S * fs + S * em .Knowing the effective strengths allows us estimate the mean pond coverage after a period of time without having to run the model.
Roughness has a different effect on freeboard sinking and enhanced melting.Freeboard sinking is roughly independent of time and proportional to the square of nondimensional roughness.Therefore, it is very sensitive to variations in roughness: doubling the ice roughness roughly quadruples the mean pond coverage due to freeboard sinking after some time.Enhanced melting depends roughly linearly on roughness.However, as roughness tends to zero, the effective strength remains nonzero, S * em ( σ → 0) → S em .Therefore, ponds on smooth ice grow primarily due to enhanced melting.Effective strength also depends on the nondimensional time, t, and is higher and more sensitive to variations in roughness early in the melt season.
Multiyear ice topography shown in Fig. 2a, dashed line, has σ ≈ 0.25 and is significantly smoother than first-year ice topography shown in Fig. 2a, solid line, which has σ ≈ 0.55.From Eq. ( 29) it follows that freeboard sinking on multiyear ice is roughly 5 times less efficient in growing the ponds than on first-year ice.

Analyzing the 0-D model yields useful insight into factors influencing the pond evolution
Extracting the dependence of a desired property on physical parameters and understanding its scaling is the main strength of our model.These types of relationships would be difficult to obtain in a more complex model.The parameters S * bi , S * mp , S * bot , and S * em control the mean rates of pond growth by melting different regions of ice.Roughly, they represent the amount of pond growth per unit time by freeboard sinking due to melting bare ice, freeboard sinking due to melting ponded ice, freeboard sinking due to melting ice bottom, and enhanced melting.Knowing these parameters allows us to estimate mean pond coverage after a period of time with significant accuracy without having to run the numerical model.Moreover, analyzing them can yield useful insight into the behavior of melt ponds under general circumstances.
We can estimate the change in magnitude of the strength of each of the growth modes when a physical parameters p changes by p as where S * i is the change in magnitude of the effective strength of the ith growth mode.This equation holds so long as the change in the physical parameter is not too large.A change in pond growth rate can then be estimated as S * = i S * i .Then, using Eq. ( 28), we can roughly estimate a change in mean pond fraction, < x >, after some time, t, following a change in physical parameter, p, as < x >≈ 1 2 S * t.This provides a means to estimate changes in mean pond coverage under different environmental conditions.
6.1 Ponds are more sensitive to changes in bare ice albedo than changes in pond albedo We will illustrate the use of effective strengths using an example where we vary the ice and pond albedos.If the bare ice albedo changes by α bi , the change in growth rate would be roughly The Cryosphere, 11, 1149-1172, 2017 If the melt pond albedo changes by α mp , the change in growth rate would be roughly It follows from these estimates that after a month the mean pond fraction would differ by roughly 4.5 % for a bare ice albedo difference of 0.1 and by around 1 % for a pond albedo difference of 0.1.Therefore, variation in pond albedo affects pond evolution roughly 5 times less than variation in bare ice albedo.This explains our observation from Fig. 6c that pond evolution is much more sensitive to variations in bare ice albedo than to variations in pond albedo.In this way, we also extract the dependence of sensitivity on physical parameters.A major difference between the two sensitivities is their dependence on the initial pond coverage: the sensitivity to pond albedo is proportional to x i , whereas the sensitivity to bare ice albedo is proportional to 1 − x i .In the above example we used x i = 0.2, which explains most of the large difference between the two sensitivities.If the pond coverage were higher, variations in the pond albedo could become more important than variations in bare ice albedo.For example, assuming no enhanced melting, the sensitivity to pond albedo would become greater than the sensitivity to bare ice albedo at 50 % pond coverage

Under global warming, pond feedback could lead to significant ice thinning
We now use the effective strengths to roughly estimate the impact of global warming on the pond coverage.At high latitudes, feedbacks due to changes in albedo, the atmospheric lapse rate, and clouds can amplify the forcing due to global warming (Holland and Bitz, 2006).For this reason forcing at high latitudes is generally larger than direct radiative forcing due to an increase in CO 2 concentration.In a global warming scenario, the pond growth rate would increase because the ice melts faster but also because ice at the beginning of the melt would be thinner.We can emulate a global warming scenario by increasing the flux F r by a certain amount, F r , and by assuming that the initial ice thickness decreases by H ≡ ∂H ∂F r F r , where ∂H ∂F r is the ice thinning per 1 W m −2 of warming.Therefore, we split the change in pond growth rate due to global warming, S * , into a contribution from direct forcing, S * F , and a contribution from ice thinning, S * H . Using the above formalism, we find The numbers in Eq. ( 33) were obtained using the default values of the parameters, and ∂H ∂F r = −0.05m 3 W −1 roughly estimated using the Eisenman and Wettlaufer (2008) model.This means that after a month's growth global warming would increase mean pond coverage by roughly 1.2 % per 1 W m −2 of warming.Nearly half of this increase in the mean pond coverage comes from an increase in the strength of enhanced melting due to ice thinning.Simulating a 30-day melt numerically using our model predicts an increase in mean pond coverage with forcing at a rate of 1.5 % per 1 W m −2 of warming for small forcing ( F r ≈ 0), which confirms the approximate validity of our linearization.For larger forcing, the sensitivity of pond coverage to forcing increases because the ice thins.Our linearized estimate, Eq. ( 33), also gives the dependence of the sensitivity on physical parameters.In a likely scenario where the forcing is around 10 W m −2 , our estimate predicts that after a month mean pond coverage would increase by around 15 %, which corresponds to around 12 cm of ice thinning solely due to the pond feedback.Ice thinning after a month directly due to forcing would be only around 9 cm, meaning that the pond feedback must be taken into account to understand ice thinning under global warming.Increased forcing could also lead to changes in initial pond coverage, changes in ice roughness, or changes in s or k.We ignored these feedbacks, as we have no way of reliably estimating ∂p ∂F r for these parameters.

Different growth modes yield different pond evolution
Each of the four growth modes has different effects on the pond coverage.We will now look in detail at each of the growth modes, their effect on the pond evolution, and their scaling with physical parameters.Figure 8 shows the dependence of growth rate on pond fraction and solutions to   26) when only one of the growth modes is active.The x axis shows the normalized time, where 0 corresponds to the beginning of the melt and 1 to entire floe being flooded.
Eq. ( 26) when only one of the strengths is nonzero, assuming a first-year ice topography.Figure 9 shows the evolution of pond coverage distribution when only one of the strengths is nonzero.
All modes of growth depend in the same way on the bulk ice density, ρ b .Each of the strengths is inversely proportional to ρ b , meaning that ponds grow faster on ice with a lower bulk density.The effect is, however, modest: within a reasonable range of 916 kg m −3 > ρ b > 750 kg m −3 , pond growth rate can vary by at most 20 %.
We will first discuss freeboard sinking.Common to all modes of freeboard sinking is the dependence on ice thickness.Each freeboard sinking growth mode is inversely proportional to the ice thickness, S * fs ∝ 1 H , meaning that, all else equal, ponds grow proportionally slower on thicker ice.
Although ice roughness may have a different effect on each of the individual modes of freeboard sinking, for simplicity we will assume that they are all affected by roughness in the same way, as parameterized in Eq. ( 29).In that case, each of these strengths is roughly proportional to the square of the nondimensional ice roughness, S * fs ∝ σ 2 , meaning that pond growth due to freeboard sinking is suppressed on smooth ice.
We will now focus on individual components of freeboard sinking.The parameter S * bi controls pond growth by freeboard sinking due to melting bare ice.On first-year ice, owing to the shape of the hypsographic curve, the pond growth rate by bare ice melting increases up to a certain pond cov-erage and decreases afterwards (Fig. 8, blue line).S * bi is proportional to the flux F bi and depends on the initial pond coverage as S * bi ∝ (1 − x i ) 2 .The quadratic dependence on initial bare ice fraction means that ponds on floes with less initial pond coverage grow faster.It also means that floes that start off less ponded can at some point become more ponded than floes that start off more heavily ponded.We can see this in Fig. 9a, where the pond coverage distribution narrows up to a certain point, after which it starts to widen again because floes with lower x i overtake the floes with higher x i .Using the default values of physical parameters of F bi = 85 W m −2 , H = 1.5 m, x i = 0.2, and σ = 0.55, we get S * bi ≈ 0.13 month −1 .The parameter S * mp controls pond growth by freeboard sinking due to melting ponded ice.The pond growth rate increases with pond fraction from 0 at x = 0 to very high values at high pond coverage and can be the dominant mode of pond growth if the pond coverage is high enough (Fig. 8, green line).For this reason, giving a representative number to pond growth rate, such as S mp , is only meaningful if the melt season is short enough such that pond coverage during that period does not change substantially.The dependence on initial pond coverage is S * mp ∝ x i (1 − x i ).For this reason the pond coverage distribution widens over time when S * mp is dominant (Fig. 9b).Using F mp = 171 W m −2 and other parameters the same as above, we get S * mp ≈ 0.07 month −1 .Although in this case, melting ponded ice affects pond evolution less than bare ice melting, it can become stronger if  26) when only one of the growth modes is active.Red curves represent the initial pond fraction distribution, blue curves represent the pond fraction distribution after a time, t, while the green curves represent the pond fraction distribution after 2t.A time used in panel (a , and in panels (c) through (f) the pond coverage is higher.For example, S * mp and S * bi are roughly the same at x = 0.35, while at x = 0.5 S * mp is roughly twice as large as S * bi .The parameter S * bot controls pond growth by freeboard sinking due to melting of the ice bottom.The pond growth rate due to bottom melting increases with increasing melt pond fraction, although more gradually than in the ponded ice melting case (Fig. 8, red line).Since the growth rate is proportional to the bare ice fraction, S * bot ∝ (1−x i ), the pond coverage distribution gets concentrated over time (Fig. 9c).Using F bot = 20 W m −2 and other parameters the same as above, we get S * bot ≈ 0.04 month −1 .The contribution from ice bottom melting becomes larger than the contribution from bare ice melting only at high x.
Now, we will turn to enhanced melting.The parameter S * em controls pond growth by enhanced melting and is the least constrained in our model due to the many poorly constrained physical processes that potentially contribute to it.Here we will only consider enhanced melting due to height-dependent processes (Eq.24) and leave lateral melting for the discussion (Sect.7.1).
Because the growth rate by enhanced melting is inversely proportional to the hypsographic curve, pond growth by enhanced melting is very fast at the beginning of the melt and decelerates afterwards (Fig. 8, cyan line).The enhanced melting strength is inversely proportional to the square of the ice thickness, S * em ∝ 1 H 2 , meaning that it is significantly more sensitive to variations in thickness than freeboard sinking.However, it is significantly less sensitive to variations in ice roughness (Eq.29).Even on perfectly smooth ice, σ = 0, ponds will grow due to enhanced melting.In that case, however, lateral melt, rather than height-dependent enhanced melting, may dominate.
The strength of enhanced melting is proportional to the height below which enhanced melting is operational, S * em ∝ s.If we take ice wetting as a physical example, this means that enhanced melting is sensitive to microphysical processes that determine how high above sea level the ice will be wet.The dependence on the parameter k depends on its magnitude.It appears in S * em in the term k−1 ds em /ds fs +1 .The term ds em /ds fs is proportional to k − 1.Therefore, if ds em /ds fs 1, enhanced melting is proportional to k − 1.However, if ds em /ds fs 1, enhanced melting becomes independent of k.Using default parameters, we find this transition happens at around k ≈ 1.2.In the example of ice wetting, this means that enhanced melting is sensitive to albedo variations near sea level when ice near sea level has a similar albedo to the rest of the floe.However, if the albedo near sea level is significantly lower than the average, pond growth is insensitive to variations in properties of ice near sea level.
Enhanced melting is proportional to the cube of the bare ice fraction, S * em ∝ (1−x i ) 3 , making it very sensitive to variations in initial pond coverage.For this reason, the pond coverage distribution gets quickly concentrated (Fig. 9d), and it is possible for initially less ponded floes to overtake initially more ponded floes.If we assume ice wetting is the only physical process responsible for enhanced melting, we can place a rough estimate on S * sm .Taking k = 1.7, s = 0.06 m, and t = 30 days, we get for default parameters S * em ≈ 0.31 month −1 .This suggests that the contribution to mean pond coverage from enhanced melting is slightly larger than the contribution from freeboard sinking after 30 days of melt.
The black line in Fig. 8 shows the total pond evolution using the default physical parameters.The pond growth rate when both freeboard sinking and enhanced melting occur is not simply a sum of the growth rates of the four modes since the equations for freeboard sinking and enhanced melting are solved separate of each other.Therefore, the dependence of growth rate on pond coverage (Fig. 8a, black line) was obtained by finding the derivative of the pond evolution curve.The pond growth rate first decreases with pond fraction, indicating that enhanced melting dominates early in the Figure 10.The red curve is the results of Skyllingstad et al. (2009).The black curve is the solution to Eq. ( 34) with F lat = K lat F mp .The pond albedo and the shortwave, longwave, sensible, and latent heat fluxes used to find F mp are the same as used in Skyllingstad et al. (2009) and K lat = 1.5.A nearly perfect agreement between the two curves suggests that a single nondimensional constant, K lat , is enough to describe pond growth by lateral melting, and the complicated physics of lateral melting are important only in determining the value of K lat .
season and then increases, indicating that freeboard sinking dominates later in the season.The pond coverage distribution using realistic parameters narrows with time (Fig. 9e).Since each growth mode affects the pond coverage distribution in a distinct way, fitting both the evolution of the mean and the standard deviation of the pond coverage distribution in observational data could add constraints on the relevant strengths.Using the above values of strengths, we find that after a month of growth bare ice melting contributes to roughly 25 % of mean pond coverage, ponded ice melting contributes to around 13 %, ice bottom melting contributes to around 7 %, and enhanced melting contributes to roughly 55 %.

Lateral melting of pond walls by pond water
In our model, we focused on vertical changes in topography, and neglected pond growth by lateral melting of pond sidewalls by pond water.We will now briefly discuss this possibility.
This type of melt was the main focus of Skyllingstad et al. (2009), who carefully calculated the lateral melt rates of pond sidewalls by pond water.The red line in Fig. 10 shows their results.The rate of change of pond fraction due to a lateral melt flux F lat is where P is the total perimeter of the ponds and A is the area of the floe.If F lat is constant and the dependence of P on pond fraction is weak, pond growth is linear, which explains the roughly linear pond coverage evolution in Skyllingstad et al. (2009).In Fig. 10, black line, we solve Eq. ( 34) assuming a lateral melt flux proportional to the ponded ice melting flux, F lat = K lat F mp , where K lat is a constant.We use the same energy fluxes used by Skyllingstad et al. (2009) and estimate P A ≈ 0.1 m −1 from the aerial photographs taken during SHEBA.A nearly perfect match is obtained with K lat = 1.5.Therefore, a single constant that relates the rate of melt of ponded ice to the rate of melt of pond walls, K lat , is enough to capture the effects of lateral melting on pond growth.This suggests that the complicated physics of lateral melting can, to a large extent, be ignored.More work would, however, be needed to determine to what degree K lat varies under different circumstances.
If we ignore the topographic variation above sea level, pond growth due to enhanced melting also becomes linear (Eq.20).Therefore, lateral melting can approximately be considered a contribution to enhanced melting, S em , although it scales differently with physical parameters than the heightdependent enhanced melting (Eq.24).It is important to note that in this model lateral melt does not depend on ice thickness, H , or on initial pond coverage, x i , although, in reality, it may depend on these to some degree.For this reason, the pond coverage distribution width does not change in time, while the mean increases linearly (Fig. 9f).
It is not simple to understand the contribution of lateral melting to pond growth when both lateral and vertical melting occur simultaneously.Each point along the pond boundary can expand either by lateral melting or by vertical melting, but not by both.This is because when a point along the pond boundary melts laterally, it creates a completely vertical slope at that point.Therefore a small vertical shift will not grow the ponds, and a large vertical shift will outgrow the lateral expansion.Therefore, if pond growth due to vertical melting is strong, the contribution from lateral melting will be small.This is consistent with observations of Polashenski et al. (2012) and Landy et al. (2014), who found that on first-year ice the contribution from lateral melting is small.However, steep topography on level multiyear ice inhibits pond expansion through vertical motion and could lead to lateral melting being the dominant mode of growth.This is consistent with our findings of a large contribution from enhanced melting to pond growth on multiyear ice during SHEBA (Fig. 5).

Effects of density variations and internal melt
So far, we have assumed that all the melt occurs on either the top or the bottom surface of the ice.However, some of the melt can happen internally, in the bulk of the ice.Internal melt occurs when trapped brine pockets with high salt content expand and dilute in order to reach a thermodynamic equilibrium with the surrounding ice.This phenomenon has been reported to occur both above and below sea level.Internal melt leads to a reduction in bulk ice density, ρ b , which in turn affects pond evolution.Accounting for internal melt correctly can be quite challenging as it requires detailed knowledge of the vertical structure of internal melt and bulk density.Nevertheless, we find that although the effects of internal melt and density variation may be significant when considered individually, if considered together they are likely small.
If internal melt is uniform throughout the vertical ice column, the only effect is a gradual reduction in ρ b over the course of the melt season, slightly increasing the pond growth rate.If, however, internal melt has a vertical structure, it will create a vertically nonuniform bulk ice density which can have more complicated effects on pond evolution.Variations in bulk density and internal melt affect pond evolution in the following ways: (1) mass transported across sea level due to rigid body movement depends on the bulk density at sea level; (2) the volume of ice removed by local melt depends on the bulk ice density at the surface; (3) freeboard height depends on average bulk densities above and below sea level; and (4) internal melt induces rigid body motion by melting mass above and below sea level, without changing the ice surface.We outline the procedure to include these effects in the pond evolution model in Appendix D. The resulting equation for pond coverage evolution has the same form as Eq. ( 26), with only the strengths modified.Here, we only qualitatively discuss our findings.Pond evolution is most sensitive to the following: 1.The difference between the internal melt rate above and below sea level, e asl −e bsl , creating a rigid body motion.
Here, e asl/bsl is the energy density used for internal melting, averaged over all ice above or below sea level.More internal melt above (below) the sea level will create an upward (downward) rigid body motion of the floe, slowing down (speeding up) pond growth.
2. The difference between the bulk ice density at the surface and the bulk ice density at sea level, ρ b (h) − ρ b (0), changing the ratio of topographic change due to local melt to rigid body motion.Using default parameters, rigid body motion is upwards, slowing down pond growth.Therefore, a lower (higher) bulk ice density at the surface relative to sea level increases (decreases) the rate of local melt relative to rigid body motion, speeding up (slowing down) pond growth.
If considered as independent processes, vertical variations in bulk ice density and internal melt can significantly alter the rate of pond growth.For example, assuming ρ b (0) = 850 kg m −3 , ρ b (h) = 750 kg m −3 , and no internal melt leads to a roughly 60 % increase in the pond growth rate.However, these processes depend on each other and have the opposite effects on pond evolution.For example, a high rate of internal melt above sea level, slowing down pond growth, will lower the bulk ice density above sea level, speeding up pond growth.
Density and internal melt can be related via a differential equation, ∂ρ b (z) ∂t = − e(z) l − ∂ρ b (z) ∂z ds rigid body dt , where z is a vertical coordinate within the ice column.Assuming vertically uniform rates of internal melt above and below sea level, an approximate long-time solution to this equation yields a vertically uniform bulk density below sea level, and a linearly decreasing bulk density above sea level.This also defines a long-time relationship between the vertical profiles of internal melt and bulk ice density, e asl −e bsl = l h ds rigid body dt (ρ b (0)− ρ b (h)).Using densities from the example in the paragraph above, and the rate of internal melt obtained in this way, leads to a roughly 10 % increase in pond growth rate, significantly less than the 60 % we found when considering only the effects of vertical density structure.
A long-time effect of vertically nonuniform internal melt and density is always a significant compensation between the two, although there may be transient effects.For this reason, we believe that including a vertical structure of density or internal melt in the simple model of pond evolution model is most likely unnecessary.

Under certain conditions, ponds can stop growing
Here, we will entertain the possibility of pond growth by vertical motion of the topography stopping entirely for a period of time.This is an example of a possible transient effect of internal melting, which, although interesting, seems unlikely.
If there is enough mass removed above sea level to induce an upward rigid body motion that is able to compensate for the effects of local melting near the sea level, points near the sea level would move upwards, ds dt > 0, and pond growth would stop.This could, for example, occur if there is strong internal melting above sea level.After a time, however, high internal melt above sea level would lower the bulk ice density at the surface, thereby increasing the rate of local melt and reinitializing pond growth.
We will use an equation for ds dt that includes the effects of vertically nonuniform internal melt and bulk ice density we derive in Appendix D, Eq. (D2).Requiring that ds dt (x) > 0 for any x, we find the condition for pond growth stopping as where ρ asl/bsl is the average bulk density above and below sea level.Using the values of internal melt and bulk densities ≈ 1, we find that in order for ponds to stop growing, k has to be less than 0.85.This is unlikely as ice near the sea level likely melts faster than ice higher up.Nevertheless, if internal melt has not had enough time to adjust densities above and below sea level, it is possible that pond growth could be stopped for a time by the action of internal melt above sea level.For example, assuming the same internal melt as in the previous example but a uniform bulk ice density (ρ b (h) = ρ b (0)), pond growth would be stopped at k = 1.In this case it is likely that growth by lateral melt would take over, as Eq. ( 35) ensures only that pond growth by vertical motions is prevented.

Conclusions
We presented a simple analytical model for melt pond evolution on permeable Arctic sea ice.The model is represented by two ordinary differential equations in which the rate of change of pond coverage depends on pond coverage.The model is governed by four parameters, S bi , S mp , S bot , and S em , that control the rate of pond growth by bare ice melting, ponded ice melting, ice bottom melting, and enhanced melting.Using this model we are able to reproduce observations well.
Our main finding is that we can estimate the mean pond coverage as a function of time without running the model by using "effective strengths": S * bi , S * mp , S * bot , and S * em .Here all the physical parameters combine in a known way which permits understanding of the behavior of pond coverage under general conditions.The most important conclusions we draw from analyzing the effective strengths are as follows: 1. Ponds grow slower on smoother ice, with freeboard sinking roughly proportional to the square of the bare ice roughness and enhanced melting increasing roughly linearly with roughness.
2. Ponds respond to both freeboard sinking and enhanced melting on first-year ice and almost entirely to enhanced melting on multiyear ice.
3. The pond growth rate is more sensitive to changes in bare sea ice albedo than changes in pond albedo unless the ice is already mostly covered in ponds.
4. Under a global warming scenario, the pond feedback could lead to ice thinning comparable to thinning due to direct forcing.
5. The dependence of ice albedo on height above sea level is likely a significant control on pond evolution.
6.The pond coverage distribution over an ensemble of floes likely narrows over time.
7. Pond evolution is insensitive to small timescale variations in the forcing.
8. If freeboard sinking is suppressed by topography, lateral melting likely plays an important role, making it a significant factor on multiyear ice.9.The complicated physics of lateral melting can be summarized by a single nondimensional constant K lat that relates the lateral melt flux to the flux used for melting the pond bottom.
10.The vertical structure of density and internal melt can likely be ignored.
As melt pond coverage is one of the key controls on summer Arctic sea ice albedo, some representation of it in GCMs is necessary for predicting the future of sea ice and its impact on global climate.With the exception of enhanced melting, our model depends only on parameters that are either available in large-scale models or that can be reasonably estimated.Therefore, if stricter constraints can be placed on the strength of enhanced melting, our model may present an accurate and computationally low-cost representation of sea level melt ponds that could be used in GCMs.

Appendix D
Here, we outline the procedure to include the effects of vertically nonuniform internal melt and bulk ice density.We assume that the bulk ice density, ρ b , and the energy density used for melting the ice internally, e, have a vertical structure, ρ b (z) and e(z), where z is positive upwards, z = 0 corresponds to sea level, and z = h corresponds to ice surface.
Mass transported across sea level depends on the bulk density at the sea level, the rate of local melting depends on the bulk ice density at the surface, and the freeboard height depends on the average densities above and below sea level, ρ asl/bsl .Internal melt above and below sea level creates a rigid body motion.This is summarized as ds loc dt (r) = −k(r) where H d is the ice draft depth defined as the volume of ice below sea level divided by the area of the ice floe, e asl/bsl is the energy density used for internal melting averaged over all ice above or below sea level, and ρ b ≡ ρ bsl −ρ asl ρ bsl is the relative difference in mean bulk density above and below sea level.
With these changes, we can find the equation for pond coverage evolution straightforwardly by repeating all of the steps from Sect. 2. We first derive the equation for the vertical motion of points near the sea level ds dt = − (k − 1) where e ≡ e bsl −e asl e bsl is the relative difference in average energy density used for internal melting below and above sea level.The two terms in square brackets correspond to enhanced melting and freeboard sinking.Then we repeat the procedure to relate Eq. (D2) to the change in pond coverage.The resulting equation has the same form as Eq. ( 26), with Here, the strength of internal melting, S int , should be included in the equation for freeboard sinking.The term ds em ds fs is given by the ratio of the two terms in square brackets in Eq. (D2).The equation for pond growth, Eq. ( 26), using the above strengths, Eq. (D3), should also be supplemented with an equation for evolution of bulk density:

Figure 1 .
Figure 1.(a) Local displacement represents the movement of a point on the ice surface as a result of ice melting at that particular point.It is a function only of local ice characteristics at that point.For both local and hydrostatic displacements the positive direction is defined as upwards.(b) Rigid body displacement represents the motion of a floe as a whole in an effort to maintain hydrostatic balance because melting removes mass above or below sea level.Melting above sea level induces an upward rigid body motion of the floe, whereas melting below sea level induces a downward motion.
express dm rigid body in terms of rigid body displacement of the floe as dm rigid body asl = ρ b A bi ds rigid body , b A bi ds rigid body ,

Figure 2 .
Figure 2. Hypsographic curves showing the percentage of the sea ice surface that is lower than a particular elevation.Pond coverage on highly permeable sea ice can be inferred from here as the intersection of sea level (horizontal blue line) with the hypsographic curve.(a) A hypsographic curve measured by Landy et al. (2014) on 25 June 2011 (solid black line) and a hypsographic curve measured during SHEBA along a 100 m long "topography profile 1" on 10 July 1998 (black dashed line).The vertical dashed lines represent the pond coverage, assuming that ice is permeable.The red line represents a fit to the part of the hypsographic curve above sea level with a tangent function, Eq. (A1).(b) Adjusted hypsographic curves for different initial pond coverage and the same ice thickness.(c) Adjusted hypsographic curves for the same initial pond coverage and different ice thickness.(d) Hypsographic curves for different shape parameters, p 1 and p 2 , defined and discussed in Appendix A, Eq. (A1).Parameter p 1 controls the amount of curvature, while p 2 controls the position of the inflection point of the tangent function.

Figure 3 .
Figure 3. Explanation of different models of pond growth.Models evolve a hypsographic curve, s(x h ), above sea level to find the pond coverage evolution.Evolution of the hypsographic curve below sea level is not relevant for pond growth and, apart from the 1-D model, is not captured well in these models.(a) Freeboard sinking shifts the entire hypsographic curve downward following a displacement of ds fs .(b) Enhanced melting acts on a constant ice fraction, δ, and there is no freeboard sinking.The hypsographic curve changes only between x h = x and x h = x + δ and remains unchanged otherwise.After a time t = s(x+δ) ds em /dt pond coverage grows by δ.The 0-D model, Eq. (26), assumes that the total pond evolution is the sum of pond evolution due to such enhanced melting and freeboard sinking (panel a).(c) The 1-D model prescribes a melt rate at each point on the hypsographic curve as a function of height above sea level, ds dt (s).(d) A simplified model that assumes both freeboard sinking and enhanced melting (Appendix B).Enhanced melting occurs only below height s.After some time, the fraction of ice affected by enhanced melting, δ, becomes constant, meaning that a constant fraction model (panel b) and a constant height model are equivalent if δ and s are related appropriately.

Figure 4 .
Figure 4. (a) A comparison between pond evolution in the 0-D model and the 1-D model.The black curve represents the 0-D model.The blue, green, and red curves represent the 1-D model for different functions k(s) shown in panel (b).These different functions were chosen such that the integral parameter < S em > (Eq.25) is the same as for the 0-D model.The yellow curve represents the 1-D model where enhanced melting acts on a constant fraction of bare ice, δ, chosen according to Eq. (22).The magenta curve represents the 1-D model with pond albedo varying with depth.There is significant agreement between all of the curves, suggesting that the simplifications made in the simple model were justified.Since including variable pond albedo does not change the pond evolution significantly, this detail can be neglected when estimating the pond coverage on permeable ice.(b) The blue, green, and red lines represent functions k(s) − 1 used to run the 1-D model.

Figure 5 .
Figure 5.A comparison between measurements of pond fraction made during SHEBA along the albedo line (red line), along a topography profile (blue dots), and our model (black line).The blue dots have been shifted downward by 0.05 to make a more obvious comparison between albedo line and topography profile trends.The black dashed line is the contribution to our model from freeboard sinking and the black dotted line is the contribution from enhanced melting.Ponds grow almost entirely due to enhanced melting as a result of the steep topography of multiyear ice.

Figure 6 .
Figure 6.Numerical solutions to Eq. (26) with parameters varied around the defaults described in the text.(a) Varying initial pond coverage.Solid lines represent solutions using full time-varying fluxes, while dashed lines represent solutions using time-averaged fluxes.The two solutions are very similar, so we subsequently use only the time-averaged fluxes.(b) Varying ice thickness.Ponds grow slower on thicker floes.(c) Varying pond and bare ice albedo.Different colors represent different bare ice albedos, and full, dotted, and dashed lines represent different pond albedos.A change in bare ice albedo has a much larger effect on pond fraction than the same change in pond albedo.(d) Varying the s and k.For k = 0.8, the ponds shrink.However, pond evolution for k < 1 is not represented well in our model, so this curve serves only as an illustration.

Figure 7 .
Figure 7. Exploring the effects of sea ice roughness.(a) Pond evolution due to pure freeboard sinking for hypsographic curves with different shape parameters p 1 and p 2 .The x axis shows nondimensional time t = t (S bi +S mp +S bot ) 1−x i .Color represents normalized roughness, σ , with blue colors corresponding to small σ and red colors corresponding to large σ .Thick red solid line represents pond evolution on the measured first-year ice hypsographic curve, and the thick red dashed line represents pond evolution on the measured multiyear ice hypsographic curve.All else equal, rougher ice has a larger pond fraction.(b) Pond evolution due to pure enhanced melting for hypsographic curves with different shapes.The x axis shows nondimensional time t = tS em 1−x i .Cartoon examples of hypsographic curves and their approximate positions along the σ axis are also shown.

Figure 8 .
Figure 8.(a) Dependence of growth rate on pond coverage for different modes of pond growth.The y axis shows the growth rate, dx dt , for each of the growth modes calculated using the default parameters and x i = 0. Pond growth rate for bare ice melting (blue line) first increases up to a certain pond coverage and then decreases.Ponded ice melting (green line) increases with pond coverage from dx dt = 0 at x = 0 to very high values at high pond coverage.The ice bottom melting rate (red line) gradually increases with pond coverage.The vertical enhanced melting rate (cyan line) decreases with pond coverage.The black line represents a realistic combination of the four growth modes and shows that pond growth is dominated by enhanced melting early in the season and by freeboard sinking late in the season.The dashed magenta line represents lateral melting estimated using parameters described in Sect.7.1.(b) Solutions to Eq. (26) when only one of the growth modes is active.The x axis shows the normalized time, where 0 corresponds to the beginning of the melt and 1 to entire floe being flooded.

Figure 9 .
Figure 9.In this figure we have evolved an ensemble of 10 5 floes with varying initial pond coverage according to Eq. (26) when only one of the growth modes is active.Red curves represent the initial pond fraction distribution, blue curves represent the pond fraction distribution after a time, t, while the green curves represent the pond fraction distribution after 2t.A time used in panel (a) is t = 1 where x i is the mean pond fraction of the initial distribution and S is an appropriate strength.We show how different growth modes have different effects on the pond fraction distribution.(a) Bare ice melting first narrows the distribution and then widens it.(b) Ponded ice melting widens the distribution.(c) Bottom ice melting narrows the distribution, while the mean of the distribution increases at an increasing rate.(d) Enhanced melting narrows the distribution, while the mean of the distribution increases at a decreasing rate.(e) Using realistic the pond distribution slowly narrows and accelerates.(f) Due to lateral melting, pond coverage distribution does not change width, and the growth is linear.

Figure C1 .
Figure C1.Determining the effective strengths, S * ≡ f ( σ, t)S.Points represent estimates of the correction f ( σ, t) for each of the curves in Fig. 7 evaluated at different times t ≡ St 1−x i.The function f ( σ, t) is evaluated as f ( σ, t) ≡ 2(< x(t) > −x i )/(St).Different colors correspond to different times with black corresponding to early in the season and magenta to late in the season.Nondimensional roughness, σ , is shown on the x axis.(a) f fs ( σ, t) evaluated for the freeboard sinking curves in Fig.7a.There is no obvious dependence on t.Freeboard sinking becomes completely suppressed as roughness tends to zero.The dashed red line represents the fit to these estimates of the form f fs ( σ, t) = a σ 2 .(b) f em ( σ, t) evaluated for the enhanced melting curves in Fig.7b.There is a clear dependence on t.Enhanced melting proceeds even as roughness tends to zero.Red dashed lines are fits to these data of the form f em ( σ, t) = 1 + c(t) σ , where c(t) ≡ 2 √ t − 3 2 .
Popović and D. Abbot:Melt pond coverage on permeable Arctic sea ice Surface elevation above sea level at point r s(x h ), s(x h ), d s Hypsographic curve, nondimensional hypsographic curve, s(x h ) = s(x h ) h , and its nondimensional derivative, d s d x h = 1−x i h ds dx h ds rigid body , ds loc (r) Change in surface elevation due to rigid body motion and due to local melting at point r ds fs , ds em , ds em /ds fs Change in surface elevation due to freeboard sinking, enhanced melting, and the magnitude of their ratio dm melt/rigid body asl/bsl Change in mass above and below sea level due to ice melting or rigid body motion x, x, 1 − x Pond fraction, normalized pond fraction x = x x i , and normalized bare ice fraction, 1 − x = Fraction of ice below an elevation given by the hypsographic curve x s Fraction of ice below s x fs (t), x em (t), x lat (t) Pond coverage evolution due to freeboard sinking, enhanced melting, and lateral melting A, A bi , A mp Areas of the floe, bare ice, and melt ponds P Total perimeter of the ponds ρ w , ρ i , ρ b Densities of salt water, pure ice, and bulk ice once all the brine has drained l Latent heat of melting H , h Initial thickness of the ice and average initial freeboard height σ , σ Bare ice roughness and nondimensional bare ice roughness, σ = σ h p 1 , p 2 Shape parameters of the hypsographic curve that control the "amount of variability" of the curve and the location of the inflection point k(r) Ratio of the melt rate at point r to the average rate of bare ice melting s Height above sea level below which there is enhanced melting δ Fraction of the ice affected by enhanced melting α bi , α mp Albedos of bare ice and melt ponds F sol , F r Solar energy flux and the sum of longwave, latent, and sensible heat fluxes F bi , F mp , F bot , F lat Fluxes of energy used for melting bare ice, ponded ice, ice bottom, and lateral melting averaged over bare ice, ponded ice, ice bottom, and the pond perimeter |F | Representative values of fluxes, e.g., their time averages K lat Constant relating the flux of energy used for melting ponded ice to the flux of energy used for lateral melting S bi , S mp , S bot , S em Strengths of bare ice melting, ponded ice melting, ice bottom melting, and enhanced melting S * bi , S * mp , S * bot , S * em Effective strengths of bare ice melting, ponded ice melting, ice bottom melting, and enhanced melting, that take into account the effects of bare ice roughness S * fs Effective strength of freeboard sinking, S * fs = S * bi + S * mp + S * www.the-cryosphere.net/11/1149/2017/The Cryosphere, 11, 1149-1172, 2017 P. em The Cryosphere, 11, 1149-1172, 2017 www.the-cryosphere.net/11/1149/2017/