Ising model for melt ponds on Arctic sea ice

Perhaps the most iconic feature of melting Arctic sea ice is the distinctive ponds that form on its surface. The geometrical patterns describing how melt water is distributed over the surface largely determine the solar reflectance and transmittance of the sea ice cover, which are key parameters in climate modeling and upper ocean ecology. In order to help develop a predictive theoretical approach to studying melting sea ice, and the resulting patterns of light and dark regions on its surface in particular, we look to the statistical mechanics of phase transitions and introduce a two-dimensional random field Ising model which accounts for only the most basic physics in the system. The ponds are identified as metastable states in the model, where the binary spin variable corresponds to the presence of melt water or ice on the sea ice surface. With the lattice spacing determined by snow topography data as the only measured parameter input into the model, energy minimization drives the system toward realistic pond configurations from an initially random state. The model captures the essential mechanism of pattern formation of Arctic melt ponds, with predictions that agree very closely with observed scaling of pond sizes and transition in pond fractal dimension.


Introduction
While snow and ice reflect most of the sunlight incident on Arctic sea ice, melt ponds absorb most of it. The ponds largely control the albedo, or solar reflectance of sea ice, as well as its transmittance [1][2][3][4][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. The ponds play a critical role in ice-albedo feedback, a key mechanism in the rapid decline of the summer Arctic ice pack [6]. In fact, by accounting for ponds in climate simulations, predicted ice pack volumes are significantly lower [7], and the yearly Arctic sea ice minimum can be accurately forecast from melt pond area in spring [4]. The impact of melt pond evolution extends into the biosphere as well [8][9][10], where the ponds act as windows for light to shine into the upper ocean, affecting Arctic marine ecology. Typical pond configurations are shown in figure 1(a).
There has been significant progress on numerical models of melt pond evolution [7,[3][4][5], although current generation melt pond parameterizations in climate models track melt water volume, not how melt water is distributed on the ice surface. However, the geometry of melt ponds and their spatial patterns impacts various sea ice and upper ocean processes such as albedo evolution, the break-up of floes, the evolution of the floe size distribution, and the patterns of light in and under the ice, which can affect photosynthetic activity and the ecology of microbial communities.
There are two key, benchmark observations of melt pond geometry that must be accounted for by a statistical physics theory of melt ponds. The first dates back to the 1998 SHEBA expedition and the measurement of melt pond sizes from images taken from helicopters [11]. The pond size distribution function prob(A) exhibits power law scaling with the observed value of the exponent ζ for pond areas in the range 10 m 2 <A<1000 m 2 being about −3/2.
Area-perimeter analysis of images of melt ponds from SHEBA as well as the 2005 Healy-Oden Trans Arctic Expedition (HOTRAX) has shown that as the ponds grow and coalesce into much larger connected structures they display a transition in fractal geometry [12], evolving from simple Euclidean shapes into complex, selfsimilar regions whose boundaries behave like space-filling curves. The fractal dimension of the boundary curves transitions from 1 to about 2 around a critical area of about 100 m 2 . In addition to constraining the geometry of melt pond evolution, the area-perimeter relationship is key to quantifying components of pond growth, such as vertical versus lateral melt, regulating the extent of the water-ice interface where lateral expansion of the ponds can occur.
Recent work shows that these geometrical characteristics are consistent with behavior exhibited by continuum percolation models [13][14][15]. In [16] a melt pond boundary is the intersection of a random surface representing the snow topography with a horizontal plane representing the water level. As the plane rises the ponds grow and coalesce. An autoregressive class of anisotropic random Fourier surfaces with correlation parameters based on snow data provides topographies that yield realistic ponds, the observed transition in fractal dimension, and a framework to analyze how the shape of the fractal transition depends on topographic characteristics.
In [17] a void model for melt ponds is introduced, where disks of varying size which represent ice and can overlap are placed randomly on the plane, with the voids between them representing the ponds. Data on pond sizes, area fractions, and correlations measured from helicopter photos of melt ponds are incorporated through parameters input into the model. The model yields the observed fractal transition and pond size distribution, and can be used to explore the generality of the behavior.
Here we address the challenge of developing a predictive theoretical model of melt ponds which accounts for the most basic physics of the system, and which yields realistic pond configurations obtained through minimization of the energy of the model. After all, we are interested in a solid-liquid phase transition from sea ice to sea water, albeit over large length and time scales. We turn then to the statistical mechanics of the Ising model to introduce such an approach [13,18]. Only the most essential physics is incorporated-in the same way that the original Ising model includes only the most basic aspects of a ferromagnet in an external magnetic field.
We envision a square lattice of surface patches or pixels of melt water or ice, corresponding to the classical spin up or spin down states, respectively. They are collectively influenced by an external forcing field, and interact only with their nearest neighbors. The energy of the melting sea ice system is expressed similarly to how the energy of a ferromagnet is estimated in the Ising Hamiltonian. Pond-like configurations, or connected regions of 'up spins,' result from a series of energy reducing updates of an initially random state. Glauber spin flip dynamics guide the flow of configurations toward realistic melt pond states which are local energy minima, or metastable states. We remark that while we can estimate the time scale associated with a spin flip-that is, melting or freezing a surface patch under certain conditions, we are not using the present model to directly describe the time evolution of ponds over the melt season.
Our introduction of a melt pond Ising model addresses a central issue in climate science, that is, linkage of scales. How can knowledge of local interactions be used to predict macroscopic behavior relevant to large scale, coarse-grained models? This is the type of fundamental problem that is solved in statistical physics [18,13] and homogenization for composite materials [19,15]. Illustrating the potentially broad applicability of this Here each site i is assigned a pre-melt ice height h i , and colored dark blue for water (s i =+1) and white for ice (s i =−1). Site P, to be updated, is adjacent to two water sites A and D, and two ice sites B and C. Since water tends to fill troughs, we require that s P =+1 if h P <0, and −1 otherwise.
approach, an Ising model for tropical convection was developed [20] to represent atmospheric processes unresolved by coarse scale climate models.

Theoretical framework
First, we recall the most general form of the classical Ising free energy where i ranges over a two dimensional square lattice with periodic boundary conditions, the s i are spin variables taking the values +1 or −1 corresponding to spin up or spin down, and i j , á ñ denotes nearest neighbors. The parameters H i and J ij represent the external magnetic field and coupling constants, respectively. In our melt pond Ising model the state variable is a binary (or spin) variable s i such that s i =+1 corresponds to absorptive melt water on the surface of our pixelated model sea ice floe and s i =−1 corresponds instead to reflective ice or snow on the surface. In addition, a temperature T can be defined which characterizes the strength of thermal fluctuations, but here we set T=0 assuming for simplicity that environmental noise does not significantly influence melt pond geometry. The two dimensional nature of the Ising model we consider here is most relevant to thinner, flatter first year sea ice, rather than thicker multiyear ice where it may be more important to include three dimensional effects.
To describe nontrivial spin clustering at zero temperature, the H i and/or J ij are chosen as random variables; the resulting models are collectively known as disordered Ising models [21]. In particular, one recovers the classical random field Ising model (RFIM) if the H i are independent random variables and the J ij are constant. At zero temperature, the system is usually assumed to follow Glauber single spin flip dynamics [22]: at each update step, the flip is accepted if  decreases and rejected if  increases. The spins are updated until no spin flip can further decrease . At this point, the system has found a local minimum of , known as a metastable state. Note that this state is not necessarily the ground state, which is the global minimum of .
Metastable states are especially relevant to physical systems near phase transitions, including supercooled liquids [23] and atmospheric aerosol particles [24]. On a short time scale, the system appears to be at an equilibrium state, but on longer time scales, it undergoes transitions between different metastable states [25]. For disordered Ising models, metastable states have been realized experimentally in, for example, doped manganites [26] and colossal magnetoresistive manganites [27]. Despite their importance, there are many unresolved issues concerning metastable states [22], with analytical results largely restricted to onedimension [28].

Random field Ising model
The key factor controlling melt pond configurations is the pre-melt sea ice topography, represented by random variables h i . 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 a=1 m is specified as the length scale above which important spatially correlated fluctuations occur in the power spectrum of sea ice topography (see supplementary materials available online at stacks.iop.org/NJP/21/063029/mmedia). We use the following update rule for Glauber dynamics, depending on whether there is a majority among the four nearest neighbors of a chosen site. If a majority exists, the site is updated to align with the majority due to 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 figure 1(b). This update rule does not depend on any parameters other than h i .
The above update rule can be restated as minimizing the classical RFIM free energy [29,30] h H s Js s, 2 with the uniformly applied field H=0 and the coupling constant J ;  +¥ see supplementary materials for a brief discussion of the H 0 ¹ case. To facilitate comparison with geophysical observations, the order parameter will be taken as the pond fraction F, which is defined as the fraction of up spins and therefore related to the magnetization M by F=(M+1)/2. At J=0, there is a unique metastable state, namely the ground state, given by s i =+1 if h i <H, and s i =−1 if h i >H. This process can only yield the correct melt pond geometry if the random variables h i are highly correlated [16]. As J increases, metastable states appear [31] at a wider range of pond fractions, with the entire range F 0, 1 Î [ ] covered for large enough J. As J  +¥, the two ground states are given by s i =+1 or s i =−1 for all i.

Geometry of metastable states
Below we present numerical results for the zero temperature Glauber dynamics of the RFIM, with the lattice size taken to be 1024×1024. The input spin configurations s i are independent binary variables (Bernoulli trials) that equal +1 with probability F in and −1 with probability F 1 in -, where F in denotes the input pond fraction. Note that these variables are uncorrelated with the h i . Following a random update sequence, the Glauber dynamics eventually yield a metastable state with output pond fraction F out . Note that this metastable state is generically distinct from the two ground states unless F in =0 or 1. Figure 2 shows the output configurations with F out =0.15, 0.30, and 0.45, which respectively result from F in =0.34, 0.42, and 0.48. This metastability is consistent with previous findings from a dynamical systems analysis [32].
The visual resemblance between the simulations in figures 2(a), (c) and the photos in figure 1(a) is now apparent, particularly in the well developed ponds where the minimum energy configurations of the model are quite evolved, coarse-grained and 'pond-like' in comparison to the purely random initial states. In the following we will analyze in detail the up spin clusters in figure 2(c) at F out =0.45. Figure 3(a) shows the log-log plot of the perimeter P versus the area A for these clusters (shown in physical units as Pa and Aa 2 ). Figure 3(b) shows the pond size distribution function prob(A). It exhibits power law scaling A A prob~z ( ) with the exponent ζ=−1.58±0.03 for pond areas in the range 10 m 2 <A<1000 m 2 , in excellent agreement with the observed value [11] of about −3/2.
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 D, defined in terms of the perimeter P and the area A as P A D . The input spin configuration corresponds to a site percolation process with occupation probability F in <0.5, below the site percolation threshold of about 0.593 [33]. The Ising model takes these purely random states as input and produces the metastable states represented by the cloud of points in figure 3(a). The upper edge of this cloud has an almost constant fractal dimension close to the theoretical value of 91/48≈1.896 for site percolation clusters right below the percolation threshold [33]. Therefore, this upper edge represents the unphysical clusters reminiscent of the original input, which are least affected by the Glauber dynamics. To identify the physical clusters that resemble real melt ponds, we thus choose the lower edge, or equivalently the smallest possible P for each A, as highlighted in figure 3(a). Within this data set, we further exclude both the smaller ponds with A<15 m 2 which are affected by the discreteness of the lattice, and the larger ponds with A>400 m 2 which are subject to substantial sampling variability because of their rareness. Figure 3(c) compares our Ising model D(A) function (thin solid black curve) with the observed fractal dimension dependence on area for real melt ponds (thick gray data curve) [12]. The model thin black curve is a best fit to the data points in the (A, P)-plane for model ponds, as in [16]. From this best fit curve we find that the transition happens around the inflection point A a 90 c 2 » m 2 . This predicted value agrees well with the observed value [12] of about 100 m 2 , with the full observed D(A) for real ponds reproduced in figure 3(d). The width of the transition regime in A log( )in figure 3(c) also agrees well with figure 3(d). Finally, supplementary figure 2 displays another quantifier of the onset of complexity that accounts for the entire collection of points in the (A, P)-plane. It yields the same critical transition area as before.

A scheme for more realistic pond boundaries
One discrepancy with observations is that our smaller model ponds are non-Euclidean on average, namely that they have an average fractal dimension greater than 1 (see figure 2(c)). To address this issue and better describe the physical process of melt pond formation, we can allow the ice topography h i to co-evolve with the spin configuration s i . A possible evolution scheme is outlined next.
Let us introduce a discrete time index n, and denote the ice topography and the spin configuration at time n respectively by h i n and s i n . The evolution from time n to time n+1 proceeds as follows.
where the function f and the random field g i n 1 + represent the deterministic and stochastic mechanisms of the topography evolution, encompassing internal processes of melting and freezing, as well as external influences such as environmental forcing, drainage processes, seasonal patterns, etc. In this evolution scheme, the system transitions between metastable states of an evolving free energy landscape, with the equilibration time estimated to be 4-5 d (see supplementary materials).
Here, instead of proposing a realistic expression for the function f, we simply consider f=0 for illustration purposes. In this case, the h i n at successive time steps n=0, 1, 2,L are independent (in both space and time) Gaussian variables with zero mean and unit variance. As shown in figure 4, the boundaries become smoother as n increases. As a result, the fractal dimensions of the smaller ponds become closer to 1, while for the larger ponds it remains close to 2, as is evident from comparing figures 3(a) and 4(e). The shapes of the simulated ponds in figure 4(b) closely resemble those of the observed melt ponds in figure 4(d). The power law scaling exponent of the pond size distribution function is found to be ζ=−1.71±0.02, as shown in figure 4(f).

Discussion
Our melt pond Ising model-with only one measured input parameter-produces ponds that are not only quite realistic in appearance, but with geometrical characteristics that quantitatively match very closely the observations on pond sizes and fractal dimension. This one parameter sets the length of a side of a square pixel in the lattice, and represents the scale above which the variations in snow topography are significant. Moreover, as  [12]. In panels (a)-(c), A and P are shown in physical units with the lattice constant a=1 m, and the number of sites is increased to 8192×8192 to improve the statistics. energy is minimized via Glauber dynamics the model creates order from disorder, flowing from a random initial state to a configuration with long range order.
The description of complex 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 climate models. Our approach based on energy minimization in statistical mechanics potentially opens new avenues for incorporating ponds, particularly in higher resolution, micro-and meso-scale models for regions up to hundreds of kilometers across. Efficient numerical algorithms which yield not only melt water volume but fast, accurate information about how it is distributed-based on the ambient conditions, would be broadly useful in sea ice dynamics, thermodynamics, and ecology. Assumptions about melt pond spatial structure influence the sub-grid scale spatial pattern of melt pond depths, meaning how water is distributed over the sea ice thickness distribution. These variations in water depth in turn markedly impact grid scale albedo.
The basic model presented here can be augmented to incorporate more detailed processes, such as the effect of changes in snow topography-potentially relevant in a changing climate. For example, effects of anisotropy in the topography can be included, as was studied in detail in [16]. The melt pond Ising model also offers the potential for efficient yet geometrically sophisticated parameterizations of melt ponds and their impact in climate models, as well as more refined models of sea ice physics and biology. In addition, the statistical physics approach developed here may be generalizable to other systems near the transition point between ice and water, such as tundra permafrost lakes, where the melting front has been described using a curve-shortening flow [34].
Minimal models such as the RFIM necessarily have limitations. Mathematically, the geometry of a fractal cannot be fully captured by its interpolation on a lattice. Physically, the RFIM is inherently unable to resolve processes at length scales smaller than the lattice constant. There, one may expect narrow water channels responsible for connecting smaller ponds into larger ponds. The inability to resolve such features likely causes the percolation threshold of the RFIM to differ from observations. For the metastable states obtained from random inputs, the percolation threshold is very close to 0.5 at H=0 (see supplementary materials). This threshold decreases as H decreases, but likely always exceeds the value for real melt ponds.