# Ising model for melt ponds on Arctic sea ice

Perhaps the most iconic feature of melting Arctic sea ice is the formation of distinctive, complex ponds on its surface during late spring. The evolution of melt ponds and their geometrical characteristics determines the albedo of sea ice, a key parameter in climate modeling [1, 2, 3, 4, 5]. However, a theoretical understanding of this evolution, and predictions of geometrical features, have remained elusive. To address this fundamental problem in polar climate science, here we introduce a two dimensional random field Ising model for melt ponds. The ponds are identified as metastable states [6, 7, 8] of the system, where the binary spin variable corresponds to the presence of melt water or ice on the sea ice surface. With only a minimal set of physical parameters, the model predictions agree very closely with observed power law scaling of the pond size distribution [9] and critical length scale where melt ponds undergo a transition in fractal geometry [10].

While snow and ice reflect most incident sunlight, melt ponds absorb most of it. The ponds largely control solar reflectance and transmittance of sea ice [1, 2, 3, 5], which in turn impact the heat and mass balances of the ice cover and the partitioning of energy in the upper ocean and lower atmosphere. Typical pond configurations are shown in Fig. 1(a). It has been found [4] that if a melt pond parameterization is included in climate model simulations, then predicted September ice volume from 1990 to 2007 is nearly 40% lower than in simulations which do not incorporate ponds, and is in much closer alignment with observations. Moreover, the yearly Arctic sea ice minimum can be accurately forecasted from melt pond area in spring [5]. The impact of melt pond evolution extends into the biosphere as well [11, 12], where the ponds act as windows for light to shine into the upper ocean, affecting Arctic marine ecology.

There has been significant progress on the development of numerical models of melt pond evolution [1, 3, 4, 5]. However, a fundamental theory of melting sea ice which accounts for observed pond characteristics has been lacking. Here we look toward statistical physics, and the Ising model in particular [13, 14], to develop such a theory. We envision surface patches or pixels of ice or melt water as collectively influenced by an external forcing field, and interacting only with their nearest neighbors.

A central issue in climate science is linkage of scales that is, how can knowledge of small scale local interactions be used to predict macroscopic behavior relevant to large scale, coarse-grained models? This is the type of question that is addressed in statistical physics [14, 13] and homogenization for composite materials [15, 16], where powerful methods of calculating macroscopic behavior from “microscopic” laws or microstructural information have been developed. Indeed, an Ising model for tropical convection was developed [17] to represent unresolved features in the atmosphere. Here we employ such methods to represent a critically important unresolved feature in the polar marine environment.

(a) | (b) |

First, we recall the most general form of the Ising free energy,

(1) |

where ranges over a two dimensional square lattice with periodic boundary conditions, and denotes nearest neighbors. In our model the state variable is a binary (or spin) variable such that corresponds to absorptive melt water on the surface of our pixelated model sea ice floe and corresponds instead to reflective ice or snow on the surface. The parameters and represent the external magnetic field and coupling constants, respectively. In addition, a temperature can be defined which controls the strength of thermal fluctuations, but here we set assuming that environmental noise does not significantly influence melt pond formation.

To describe nontrivial spin clustering at zero temperature, the and/or are chosen as random variables; the resulting models are collectively known as disordered Ising models [18]. In particular, one recovers the classical random field Ising model (RFIM) if the are independent random variables and the are constant. At zero temperature, the system is usually assumed to follow Glauber single spin-flip dynamics [19]: at each update step, the flip is accepted if decreases and rejected if increases. The system eventually converges to a local minimum of , known as a metastable state.

Metastable states are especially relevant to physical systems near phase transitions, including supercooled liquids [20] and atmospheric aerosol particles [21]. For disordered Ising models they have been realized experimentally in, for example, doped manganites [22] and colossal magnetoresistive manganites [23]. Despite their importance, metastable states are not completely understood theoretically [19], with analytical results largely restricted to 1D [24] and many intricate issues remaining in 2D [25].

The key factor controlling melt pond configurations is the pre-melt ice topography, represented by random variables . In the spirit of creating order from disorder, these variables are assumed to be independent Gaussian with zero mean and unit variance. The lattice constant m is specified as the length scale above which important spatially correlated fluctuations occur in the power spectrum of sea ice topography (see Supplementary Methods). We use the following update rule for Glauber dynamics, depending on whether there is a majority among the four neighbors of a chosen site. If a majority exists, the site is updated to align with the majority because of heat diffusion between neighboring sites. Otherwise, we introduce a tiebreaker rule that describes the tendency for water to fill troughs: the chosen site is updated to ice if its pre-melt ice height is positive, and water otherwise; see Fig. 1(b). Note that this update rule does not depend on any parameters other than .

The above update rule can be restated as minimizing the classical RFIM free energy [6, 7, 8]

(2) |

with the uniformly applied field and the coupling constant ; see Supplementary Methods for a brief discussion of the case. To facilitate comparison with geophysical observations, the order parameter will be taken as the pond fraction , which is defined as the fraction of up-spins and therefore related to the magnetization by . At , there is a unique metastable state given by if , and if . This process can only yield the correct melt pond geometry if the random variables are highly correlated [26]. As increases, metastable states appear [27] at a wider range of pond fractions, with the entire range covered for large enough .

Below we present numerical results for the zero temperature Glauber dynamics of the RFIM, with Monte Carlo steps used for each simulation and the lattice size taken to be . The input spin configurations are independent binary variables that equal with probability and with probability , where denotes the input pond fraction. Note that these variables are uncorrelated with the . Following a random update sequence, Glauber dynamics eventually yield a metastable state with output pond fraction . Fig. 2 shows the output configurations with , , and . This metastability is consistent with previous findings from a dynamical systems analysis [28].

(a) | (b) | (c) |

The up-spin clusters in Fig. 2(c) at correspond to well developed melt ponds [10]. Fig. 3(a) shows the log-log plot of the perimeter versus the area for these clusters (shown in physical units as and ). Fig. 3(b) shows the pond size distribution function . It exhibits power law scaling with the exponent for pond areas in the range m m, in excellent agreement with the observed value [9] of about .

(a) | (b) | (c) | (d) |

A key feature of multi-cluster systems is the tendency for smaller clusters to have simple shapes and larger clusters to have complex shapes. This onset of complexity can be quantified by an increase in the fractal dimension , defined in terms of the perimeter and the area as . To find the critical area above which shapes do not remain simple, we choose the smallest possible for each , or equivalently the lower edge of the cluster of points in the -plane. Fig. 3(c) shows the function computed for our model, illustrating the fractal dimension transition from to around a critical area . Fitting a suitable smooth function to the data points [26], we find that the transition happens around the inflection point m. This predicted value agrees well with the observed value [10] of about 100 m, as reproduced in Fig. 3(d). The width of the transition regime in in Fig. 3(c) also agrees well with Fig. 3(d). Finally, Supplementary Fig. 2 displays another quantifier of the onset of complexity that accounts for the entire cluster of points in the -plane. It yields the same critical area as before, m.

Minimal models such as the RFIM necessarily have limitations. In particular, the RFIM has a percolation threshold very close to 0.5 at (see Supplementary Methods). This threshold decreases as decreases, but likely always exceeds the value for real melt ponds. This discrepancy may be attributed to unresolved processes at smaller scales, and/or the observed pre-melt ice topography being spatially correlated rather than completely random (see Supplementary Fig. 1). We anticipate that, based on a significant amount of observational data, a detailed scheme for choosing the initial spin configuration and update sequence may be formulated. See also Supplementary Methods for possible modifications to the update rule with an alternative free energy yielding similar predictions.

The interpretation of complex Arctic melt ponds in terms of a simple disordered system may well advance our ability to model the future trajectory of the Arctic sea ice pack, e.g., through parameterizations in global climate models [29]. In addition, the statistical physics approach developed here may be generalizable to other systems near the transition point between ice and water, such as permafrost tundra lakes [30].

## Acknowledgments

We gratefully acknowledge support from the Division of Mathematical Sciences and the Division of Polar Programs at the U.S. National Science Foundation (NSF) through Grants DMS-1009704, ARC-0934721, DMS-0940249, DMS-1413454, and DMS-0940262. We are also grateful for support from the Office of Naval Research (ONR) through Grant N00014-13-10291. Y.M. acknowledges support from a Vice Chancellor’s Research Fellowship at Northumbria University. I.S. acknowledges support from the RFBR under the Grant #16-31-60070 mol_a_dk. Finally, we would like to thank the NSF Math Climate Research Network (MCRN) for their support of this work.

Y.M., I.S., C.S., and K.G. proposed the model. Y.M. and C.S. performed the numerical work. All authors contributed significantly to writing the manuscript.

## Competing financial interests

The authors declare no competing financial interests.

## Supplementary Methods

Lattice constant. The lattice constant must be small relative to the 10-20 m length scales prominent in sea ice and snow topography[32]. We set m as the length above which the power spectral density (psd) of observed snow topography exceeds a null red noise spectrum (Supplementary Fig. 1). For this calculation, we used 13 radar transects collected during the Surface Heat Budget of the Arctic Ocean (SHEBA) project [33]. To estimate the psd via the Welch modified periodogram, we calculated the power spectrum for each transect with a Hanning window and 50% segment overlap, and then averaged the results across the transects. We calculated the corresponding null red noise spectrum based on lag-one spatial autocorrelation[34] averaged across the transects.

Alternative quantifier of the onset of complexity. To account for the entire cluster of points in the -plane in Fig. 3(a), we define a new quantifier of the onset of complexity as the variance of , hereafter referred to as the elasticity. As shown in Supplementary Fig. 2, there exists a critical area such that increases with for simple ponds with , and decreases with for complex ponds with . The onset of complexity may then be identified with maximum elasticity, which occurs at m. This coincides with the critical area determined from Fig. 3(c) by the inflection point in the best fit.

Percolation threshold and correlation length exponent. For a two dimensional square lattice with occupation probability , the site-site correlation function gives the probability that a site at is a member of the same cluster as a site at . The function is assumed to decay with large distance according to

(S1) |

where is referred to as the correlation length. Theory indicates that should obey

(S2) |

where is the universal critical exponent in two dimensions and is the percolation threshold. For the two-dimensional square site lattice, [35]. For the RFIM, analysis of 5,000 model realizations on lattices yields a value close to (Fig. 3a), with correlation lengths aligning reasonably with the universal exponent (Fig. 3b). This result indicates that the spatial correlation structure of melt ponds in this model is sufficiently short-ranged so that the system falls within a standard universality class [36].

Time scale. The time scale for melt pond formation can be generally identified with the typical time taken to flip a spin in Glauber dynamics. After the RFIM decides that a spin flip is energetically favorable, we assume that the actual spin-flip process is facilitated by radiation balance [37]. The incoming shortwave radiation is , where is the mean solar insolation during polar summer, and is the surface albedo, for water and for ice. The outgoing longwave radiation is , where now is the Stefan-Boltzmann constant, and is the surface temperature, approximately K for both water and ice. Therefore the rate of heat loss for ice is , and the rate of heat gain for water is . On the other hand, the energy per unit area required for freezing a water column or melting an ice column is , where is the latent heat of fusion, is the density of water or ice (taken to be the same for simplicity), and is a realistic value for the height of the active layer. Therefore, the time intervals needed to freeze a water site or to melt an ice site are, respectively, days and days, both of which are reasonable, given this rough estimation.

Nonzero uniformly applied field. Let us choose and keep in the RFIM given by Eq. (2). Then the tiebreaker rule for a chosen site changes to if , and if , which favors ice for and water for . Here we only consider two limiting cases when the tiebreaker rule completely favors ice or water: (I) ; (II) . In these cases, the random field does not affect the kinetics, so the RFIM reduces to the classical Ising model without disorder,

(S3) |

The corresponding metastable states are known as Wulff droplets [38]. In case (I) the up-spin clusters are more elongated, and the percolation threshold is below 0.5. In case (II) the up-spin clusters are more circular, and the percolation threshold is above 0.5. These geometrical features afforded by varying (and possibly also ) provide additional prospects to describe detailed shapes of real melt pond patterns.

Alternative update rule and free energy. Let us retain the RFIM update rule when a majority exists among the neighboring sites, but adopt the following tiebreaker rule: the chosen site is updated to ice if its pre-melt ice height is larger than the average between the two neighboring ice sites, and water otherwise. For example, in Fig. 1(b) we require that if , and otherwise. This new update rule can be restated as minimizing an interfacial energy between water and ice: if a water site neighbors an ice site , then a penalty is imposed, where is a constant. The total free energy can then be written in two equivalent forms,

(S4) |

where and represent, respectively, the discrete Laplacian at site and the average between sites ,

(S5) |

The new “effective” random fields , being the curvature of , are more correlated than the by themselves. As a result, at output pond fraction , the critical area for the transition in fractal dimension and the critical area for maximum elasticity are both m. The corresponding power law scaling exponent for the pond size distribution is . It may be interesting to compare these geometrical characteristics with classical ferromagnetic random bond Ising models [39].

## References

### References

- Scott, F. & Feltham, D. L. A model of the three-dimensional evolution of Arctic melt ponds on first-year and multiyear sea ice. J. Geophys. Res. 115, C12064 (2010).
- Perovich, D. K. & Polashenski, C. Albedo evolution of seasonal Arctic sea ice. Geophys. Res. Lett. 39, L08501 (2012).
- Polashenski, C., Perovich, D. & Courville, Z. The mechanisms of sea ice melt pond formation and evolution. J. Geophys. Res. 117, C01001 (2012).
- Flocco, D., Schroeder, D., Feltham, D. L. & Hunke, E. C. Impact of melt ponds on Arctic sea ice simulations from 1990 to 2007. J. Geophys. Res. 117, C09032 (2012).
- Schröder, D., Feltham, D. L., Flocco, D. & Tsamados, M. September Arctic sea-ice minimum predicted by spring melt-pond fraction. Nat. Clim. Change 4, 353–357 (2014).
- Andelman, D. & Joanny, J.-F. Metastability in the random-field Ising model. Phys. Rev. B 32, 4818–4821 (1985).
- Grant, M. & Gunton, J. D. Metastable states in the random-field Ising model. Phys. Rev. B 35, 4922–4928 (1987).
- Perez-Reche, F. J., Rosinberg, M. L. & Tarjus, G. Numerical approach to metastable states in the zero-temperature random-field Ising model. Phys. Rev. B 77, 064422 (2008).
- Perovich, D. K., Tucker, W. B. & Ligett, K. Aerial observations of the evolution of ice surface conditions during summer. J. Geophys. Res. 107, C000449 (2002).
- Hohenegger, C., Alali, B., Steffen, K. R., Perovich, D. K. & Golden, K. M. Transition in the fractal geometry of Arctic melt ponds. The Cryosphere 6, 1157–1162 (2012).
- Arrigo, K. R. et al. Massive phytoplankton blooms under Arctic sea ice. Science 336, 1408–1408 (2012).
- Nicolaus, M., Katlein, C., Maslanik, J. & Hendricks, S. Changes in Arctic sea ice result in increasing light transmittance and absorption. Geophys. Res. Lett. 39, L24501 (2012).
- Yeomans, J. M. Statistical Mechanics of Phase Transitions (Clarendon Press, 1992).
- Christensen, K. & Moloney, N. R. Complexity and Criticality (Imperial College Press, 2005).
- Milton, G. W. The Theory of Composites. Cambridge Monographs on Applied and Computational Mathematics, Vol. 6 (Cambridge University Press, 2002).
- Torquato, S. Random Heterogeneous Materials: Microstructure and Macroscopic Properties (Springer-Verlag, 2002).
- Khouider, B., Majda, A. J. & Katsoulakis, M. A. Coarse-grained stochastic models for tropical convection and climate. Proc. Natl Acad. Sci. USA 100, 11941–11946 (2003).
- Young, A. P. Spin Glasses and Random Fields (World Scientific, 1998).
- Krapivsky, P. L., Redner, S. & Ben-Naim, E. A Kinetic View of Statistical Physics (Cambridge University Press, 2010).
- Büchner, S. & Heuer, A. Metastable states as a key to the dynamics of supercooled liquids. Phys. Rev. Lett. 84, 2168–2171 (2000).
- Rood, M., Shaw, M., Larson, T. & Covert, D. Ubiquitous nature of ambient metastable aerosol. Nature 337, 537–539 (1989).
- Moreo, A., Mayr, M., Feiguin, A., Yunoki, S. & Dagotto, E. Giant cluster coexistence in doped manganites and other compounds. Phys. Rev. Lett. 84, 5568–5571 (2000).
- Wu, W. et al. Magnetic imaging of a supercooling glass transition in a weakly disordered ferromagnet. Nat. Mater. 5, 881–886 (2006).
- Derrida, B. & Gardner, E. Metastable states of a spin glass chain at 0 temperature. J. de Phys. 47, 959–965 (1986).
- Newman, C. & Stein, D. Metastable states in spin glasses and disordered ferromagnets. Phys. Rev. E 60, 5244–5260 (1999).
- Bowen, B., Strong, C. & Golden, K. M. Modeling the fractal geometry of Arctic melt ponds using the level sets of random surfaces. J.Fract. Geom., in press (2017).
- Wu, Y. & Machta, J. Ground states and thermal states of the random field Ising model. Phys. Rev. Lett. 95, 137208 (2005).
- Sudakov, I., Vakulenko, S. A. & Golden, K. M. Arctic melt ponds and bifurcations in the climate system. Commun. Nonlinear Sci. Numer. Simul. 22, 70 – 81 (2015).
- Flocco, D., Feltham, D. L. & Turner, A. K. Incorporation of a physically based melt pond scheme into the sea ice component of a climate model. J. Geophys. Res. 115, C08012 (2010).
- Sudakov, I. & Vakulenko, S. A. A mathematical model for a positive permafrost carbon-climate feedback. IMA J. Appl. Math. 80, 811–824 (2015).
- Petrich, C. et al. Snow dunes: A controlling factor of melt pond distribution on Arctic sea ice. J. Geophys. Res. 117, C09029 (2012).
- Sturm, M., Holmgren, J. & Perovich, D. K. Winter snow cover on the sea ice of the Arctic Ocean at the Surface Heat Budget of the Arctic Ocean (SHEBA): Temporal evolution and spatial variability. J. Geophys. Res. 107, C000400 (2002).
- Gilman, D. L., Fuglister, F. J. & Mitchell, J. M. On the power spectrum of “red noise”. J. Atmos. Sci. 20, 182–184 (1963).
- Newman, M. E. J. & Ziff, R. M. Efficient Monte Carlo algorithm and high-precision results for percolation. Phys. Rev. Lett. 85, 4104–4107 (2000).
- Isichenko, M. B. Percolation, statistical topography, and transport in random media. Rev. Mod. Phys. 64, 961–1043 (1992).
- Pierrehumbert, R. T. Principles of Planetary Climate (Cambridge University Press, 2010).
- Schonmann, R. H. & Shlosman, S. B. Wulff droplets and the metastable relaxation of kinetic Ising models. Commun. Math. Phys. 194, 389–462 (1998).
- Jacobs, A. & Coram, C. Ferromagnetic random-bond Ising model: Metastable states and complexity of the energy surface. Phys. Rev. B 36, 3844–3850 (1987).