By Kenneth M. Golden, Yiping Ma, Courtenay Strong, and Ivan Sudakov
When the snow on top of Arctic sea ice begins to melt in late spring, small pools of water form on the surface. As the melt season progresses, these simply shaped meter-scale pools grow and coalesce into kilometer-scale labyrinths of cerulean blue with complex, self-similar boundaries. The fractal dimension of these boundaries transitions from one to roughly two as the area increases through a critical regime that is centered around 100 square meters [4]. While the white, snowy surface of the sea ice reflects most of the incident sunlight, the darker melt ponds act like windows and allow significant light to penetrate the ice and seawater underneath. Melt ponds thus help control the amount of solar energy that the ice pack and upper ocean absorb, strongly influencing ice melting rates and the ecology of the polar marine environment. They largely determine sea ice albedo—the ratio of reflected to incident sunlight—which is a key parameter in climate modeling.
When viewed from a helicopter, the beautiful patterns of dark and light on the surface of melting sea ice are reminiscent of structures that applied mathematicians sometimes see when studying phase transitions and coarsening processes in materials science. They also resemble the complex regions of aligned spins, or magnetic domains, that are visible in magnetic materials. Figure 1 compares two examples of magnetic domains with similar patterns that are formed by melt ponds on Arctic sea ice. Magnetic energy is lowered when nearby spins align with each other, which produces the domains. At higher temperatures, thermal fluctuations dominate the tendency of the domains’ magnetic moments to also align, with no net magnetization \(M\) of the material unless one applies an external magnetic field \(H\) to induce alignment. However, the tendency for overall alignment takes over at temperatures below the Curie point \(T_c\), and the material remains magnetized even as the applied field \(H\) vanishes, where the remaining non-zero magnetization \((M \neq 0)\) is called spontaneous or residual.
Figure 1. Comparison of magnetic domains and the patterns of meltwater on Arctic sea ice. 1a. Magnetic domains in cobalt, roughly 20 microns across. 1b. Arctic melt pond, roughly 100 meters across. 1c. Magneto-optic Kerr effect microscope image of maze-like domain structures in thin films of cobalt-iron-boron, roughly 150 microns across. 1d. Similarly-structured melt ponds, roughly 70 meters across. Figure 1a courtesy of [9], 1b and 1d courtesy of Donald Perovich, 1c courtesy of [10].
The prototypical model of a magnetic material based on a lattice of interacting binary spins is the Ising model, which was proposed in 1920 by Ernst Ising’s Ph.D. advisor Wilhelm Lenz. This model incorporates only the most basic physics of magnetic materials and operates on the principle that natural systems tend toward minimum energy states.
Consider a finite box \(\Lambda \subset \mathbb{Z}^2\) that contains \(N\) sites. At each site, a spin variable \(s_i\) can take the values \(+1\) or \(-1\) (see Figure 2). To illustrate our melt pond Ising model, we formulate the problem of finding the magnetization \(M(T,H)\)—or order parameter—of an Ising ferromagnet at temperature \(T\) in field \(H\). The Hamiltonian \(\mathscr{H}\) with ferromagnetic interaction \(J \geq 0\) between nearest neighbor pairs is given by
\[\mathscr{H}_{\omega}=-H\sum_i{s_i} - J\sum_{<i,j>}{s_is_j}\]
for any configuration \(\omega \in \Omega = \{-1,1\}^N\) of the spin variables. The canonical partition function \(Z_N\), which yields the system’s observables, is given by
\[Z_N(T,H)=\sum_{\omega\in\Omega}\exp(-\beta {\mathscr{H}}_\omega) = \exp(-\beta N f_N),\]
where \(\beta=1/kT\), \(k\) is Boltzmann’s constant, \(\exp(-\beta {\mathscr{H}}_\omega)\) is the Gibbs factor, and \(f_N\) is the free energy per site: \(f_N(T,H)=(-1/\beta N)\log{Z_N(T,H)}\).
Figure 2. Lattice models in statistical mechanics. 2a. Two-dimensional (2D) Ising model, with spins either up or down at each lattice site. 2b. Spin configuration. Spin-up sites are blue and spin-down sites are white. Image courtesy of Kenneth Golden.
The magnetization \(M(T,H) = \lim\nolimits_{N\rightarrow\infty} \frac{1}{N}\sum\nolimits_{i}{s_i}\) is averaged over \(\omega \in \Omega\) with Gibbs’ weights and expressed in terms of the free energy \(f(T,H) = \lim\nolimits_{N\rightarrow\infty} f_N(T,H)\):
\[M(T,H)=-\frac{\partial f}{\partial H}.\]
The model’s rich behavior is exemplified in the existence of a critical temperature \(T_c\)—the Curie point—where \(M(T) = \lim\nolimits_{H \rightarrow 0}M(T,H)>0\) for \(T<T_c\) and \(M(T)=0\) for \(T \geq T_c\). Universal power law asymptotics for \(M(T) \rightarrow 0\) as \(T \rightarrow T_c^-\) are independent of the lattice type and other local details.
The Metropolis algorithm is a common method for numerically constructing equilibrium states of the Ising ferromagnet. In this approach, a randomly-chosen spin either flips or does not flip based on which action lowers or raises the energy. \(\Delta E\) represents the change in magnetostatic energy from a potential flip (as measured by \(\mathscr{H}_{\omega}\)), and the spin is flipped if \(\Delta E \leq 0\). If \(\Delta E > 0\), he probability of the spin flipping is given by the Gibbs factor for \(\Delta E\). Sweeping through the whole lattice and iterating the process many times attains a local minimum in the system’s energy.
We have adapted the classical Ising model to study and explain the observed geometry of melt pond configurations and capture the fundamental physical mechanism of pattern formation in melt ponds on Arctic sea ice [5]. While previous studies have developed important and instructive numerical models of melt pond evolution [2, 3], these models were somewhat detailed and did not focus on the way in which meltwater is distributed over the sea ice surface. Our new model is simplistic and accounts for only the system’s most basic physics. In fact, the only measured parameter is the one-meter lattice spacing, which is determined by snow topography data.
The simulated ponds are metastable equilibria of our melt pond Ising model. They have geometrical characteristics that agree very closely with observed scaling of pond sizes [6] and the transition in pond fractal dimension [4]. Researchers have also developed continuum percolation models that reproduce these geometrical features [1, 8].
We aim to use our Ising model to introduce a predictive capability to cryosphere modeling based on ideas of statistical mechanics and energy minimization, utilizing just the essential physics of the system. The model consists of a two-dimensional lattice of \(N\) square patches, or pixels, of meltwater \((s_i=+1)\) or ice \((s_i=-1)\), which correspond to the spin-up or spin-down states in the Ising ferromagnet. Configurations \(\omega \in \Omega = \{-1,1\}^N\) of the spin field \(s_i\) represent the distribution of meltwater on the sea ice surface. Each patch interacts only with its nearest neighbors and is influenced by a forcing field. However, sea ice surface topography—which can vary from site to site and influence whether a patch comprises water or ice—plays the role of the applied field in our melt pond Ising model. Our model is then actually a random field Ising model, and one can write the Hamiltonian as
\[\mathscr{H}_{\omega}=-\sum_i{(H-h_i) s_i} - J\sum_{<i,j>}{s_is_j}.\]
Here, \(h_i\) are the surface heights (taken to be independent Gaussian variables with mean zero) and \(H\) is a reference height (taken to be zero in the model’s simplest form). The spin field \(s_i\) is reorganized to lower the free energy, and the order parameter is the pond area fraction \(F = (M + 1)/2\), which is directly related to sea ice albedo. We set temperature \(T = 0\) and assume for simplicity that environmental noise does not significantly influence melt pond geometry.
Figure 3. Comparison of real Arctic melt ponds with metastable equilibria in our melt pond Ising model. 3a. Ising model simulation. 3b. Real melt pond photo. Figure 3a courtesy of Yiping Ma, 3b courtesy of Donald Perovich.
Independent flips of a weighted coin determine the system’s initial random configuration. A pixel or site has a probability \(p\) of its spin being \(+1\), or meltwater. The system then updates based on simple rules: pick a random site \(i\) and update \(s_i\) as follows. If a majority exists among \(s_i\)'s four nearest neighbors, we assume that heat diffusion drives \(s_i\) to agree with this majority. Otherwise we assume water’s tendency to fill troughs, as determined by the local value of the random field \(h_i\). This update step, which corresponds to energy minimization via Glauber spin flip dynamics, iterates until \(s_i\) becomes steady. The spin-up or meltwater clusters in the final configurations of the spin field \(s_i\) exhibit geometric characteristics that agree surprisingly well with observations of Arctic melt ponds (see Figure 3). The final configuration is a metastable state — a local minimum of \(\mathscr{H}_{\omega}\). As neighboring sites exchange heat, spins tend to align to minimize energy. In doing so, they coarsen away from the purely random initial state. The emergence of this order from disorder is a central theme in statistical physics and an attractive feature of our approach.
The ability to efficiently generate realistic pond spatial patterns may enable advances in how researchers account for melt ponds and many related physical and biological processes in global climate models (GCMs). Typical GCM grid spacing is tens to hundreds of kilometers, so melt ponds are subgrid-scale and thus too small to resolve on the model grid. Instead, GCMs use parameterizations to specify a pond fraction. Specifically, modern parameterizations in GCMs track a thermodynamically-driven meltwater volume and distribute it over the sea ice thickness classes that are present in a grid cell, beginning with the thinnest class since it presumably has the lowest ice height [3]. This yields a pond fraction \(F\) and a first-order approximation to sea ice albedo, \(\alpha_{sea \; ice} = F \alpha_{water} + (1-F) \alpha_{snow}\), but does not address how the pond’s area is organized spatially. Our simple model provides a framework for prescribing a subgrid-scale spatial organization whose realistic fractal dimension or area-perimeter relation could have important influences on pond evolution [7].
At this stage, total agreement between this simple model and the real world is too much to ask. The Ising model is unable to resolve features that are smaller than the lattice constant, and the metastable state also inherits certain unrealistic features from the purely random initial condition. Nonetheless, the model may be able to use more sophisticated rules to reproduce actual melt pond evolution. We anticipate that emerging techniques—like machine learning—will deduce such evolutionary rules from observational data.
References
[1] Bowen, B., Strong, C., & Golden, K.M. (2018). Modeling the fractal geometry of Arctic melt ponds using the level sets of random surfaces. J. Fract. Geom., 5, 121-142.
[2] Flocco, D., & Feltham, D.L. (2007). A continuum model of melt pond evolution on Arctic sea ice. J. Geophys. Res., 112, C08016.
[3] Flocco, D., Feltham, D.L., & Turner, A.K. (2010). Incorporation of a physically based melt pond scheme into the sea ice component of a climate model. J. Geophys. Res., 115(C8), C08012.
[4] Hohenegger, C., Alali, B., Steffen, K.R., Perovich, D.K., & Golden, K.M. (2012). Transition in the fractal geometry of Arctic melt ponds. The Cryosphere, 6(5), 1157-1162.
[5] Ma, Y.P., Sudakov, I., Strong, C., & Golden, K.M. (2019). Ising model for melt ponds on Arctic sea ice. New J. Phys., 21, 063029.
[6] Perovich, D.K., Tucker, W.B., & Ligett K. (2002). Aerial observations of the evolution of ice surface conditions during summer. J. Geophys. Res., 107(C10), C000449.
[7] Polashenski, C., Perovich, D., & Courville, Z. (2012). The mechanisms of sea ice melt pond formation and evolution. J. Geophys. Res., 117(C1), C01001.
[8] Popović, P., Cael, B.B., Silber, M., & Abbot, D.S. (2018). Simple rules govern the patterns of Arctic sea ice melt ponds. Phys. Rev. Lett., 120(14), 148701.
[9] Portland Waldorf School. (2008). Electromagnetism. Retrieved from http://www.alchemical.org/em/index.html.
[10] Yamanouchi, M., Jander, A., Dhagat, P., Ikeda, S., Matsukura, F., & Ohno, H. (2011). Domain structure in CoFeB thin films with perpendicular magnetic anisotropy. IEEE Magn. Lett., 2, 3000304.
Kenneth M. Golden is a distinguished professor of mathematics and an adjunct professor of biomedical engineering at the University of Utah. His research is focused on developing mathematics of composite materials and statistical physics to model sea ice structures and processes. Yiping Ma is a Vice Chancellor’s Fellow of Mathematics at Northumbria University. His research focuses on nonlinear dynamics and statistical physics with applications to diverse areas, including materials, optics, climate, and fluids. Courtenay Strong is an associate professor of atmospheric sciences at the University of Utah. A substantial component of his research focuses on modeling and analysis of the cryosphere, which includes sea ice and snow. Ivan Sudakov is an assistant professor in the Department of Physics at the University of Dayton and a Kavli Institute for Theoretical Physics Scholar. He specializes in data analysis and mathematical modeling for physical and living systems.