Statistical Mechanics and the Climatology of the Arctic Sea Ice Thickness Distribution.

We study the seasonal changes in the thickness distribution of Arctic sea ice, g(h), under climate forcing. Our analytical and numerical approach is based on a Fokker-Planck equation for g(h) (Toppaladoddi and Wettlaufer in Phys Rev Lett 115(14):148501, 2015), in which the thermodynamic growth rates are determined using observed climatology. In particular, the Fokker-Planck equation is coupled to the observationally consistent thermodynamic model of Eisenman and Wettlaufer (Proc Natl Acad Sci USA 106:28-32, 2009). We find that due to the combined effects of thermodynamics and mechanics, g(h) spreads during winter and contracts during summer. This behavior is in agreement with recent satellite observations from CryoSat-2 (Kwok and Cunningham in Philos Trans R Soc A 373(2045):20140157, 2015). Because g(h) is a probability density function, we quantify all of the key moments (e.g., mean thickness, fraction of thin/thick ice, mean albedo, relaxation time scales) as greenhouse-gas radiative forcing, Δ F 0 , increases. The mean ice thickness decays exponentially with Δ F 0 , but much slower than do solely thermodynamic models. This exhibits the crucial role that ice mechanics plays in maintaining the ice cover, by redistributing thin ice to thick ice-far more rapidly than can thermal growth alone.

Introduction

Arctic sea ice is one of the most sensitive components of the Earth’s climate system and serves as a bellwether for global scale climate change. The recent decline in both the areal extent and the average thickness of sea ice, as evidenced by satellite and submarine measurements, drives study of its origins [8]. The key quantity of interest in the geophysical-scale description of sea ice is its volume; while daily areal extent is routinely measured using satellites, it is a challenge to understand the evolution of the ice volume because of the difficulties involved in the measurement of the thickness, h [8].

To study the evolution of the ice volume, one could treat ice as a continuum and construct the mass, momentum, and energy balance equations [18]. However, such a description is incomplete without the knowledge of the rheology and physical properties, such as the albedo and thermal growth rate, of the ice pack. These physical properties depend strongly on the thickness. Thus, this implies that in order to complete a continuum description, one should first determine these properties for the ice pack.

The key step in the construction of such a description was taken in 1975 by Thorndike et al. [14], who introduced the concept of thickness distribution, g(h). It is defined as follows: Consider a region with area R that is sufficiently large to contain a range ice of different thicknesses. Then the integralgives the fraction of that area (A / R) that contains ice of thicknesses between and , and the dependence of g(h) on space and time is implicit. The spatio-temporal evolution of g(h), subject to wind, thermal and mechanical forcing, is governed by [14]:where is the horizontal velocity of ice pack, f is the thermal growth/melt rate of ice, and is the redistribution function that accounts for all the mechanical interactions between ice floes (ridging, rafting, and formation of open water). The principal difficulty in solving Eq. 2 came from whose general form could not be deduced from observations. Thorndike et al. [14] separately considered the cases of the formation of open water and pressure ridges, and constructed simple models of for these events based on physical arguments. The general form of was taken to be the combination of the above mentioned cases. Numerical integration of Eq. 2 using initial conditions from the limited submarine measurements resulted in g(h)’s that were qualitatively similar to those from observations. However, Eq. 2 remained intractable due to the lack of a closed mathematical form for . Indeed, Thorndike et al. [14] noted that “The present theory suffers from a burdensome and arbitrary redistribution function .”

Thorndike, in a later study [15], made two calculations in order to understand the nature of and its role in the evolution of g(h). In the first calculation he obtained by assuming a steady state and solving for:where is the annually averaged thermal growth rate from the one-dimensional thermodynamic model of Maykut & Untersteiner (MU71) [10], and g(h) was taken from observations. Depending on the values of , the solutions displayed the following features: (a) provided a source of open water; (b) ice of thickness less than a certain value was used to build pressure ridges, and hence was a sink for this range of thickness; and (c) was a source of ice thicker than .

For his second calculation, Thorndike formulated the original equation as a Markov process; and by assuming the forms of and he solved for the steady state. In constructing the matrices of and he used the following principle: If ice of initial thickness grew either by thermal growth or by ridging to a final thickness , then the process that led to this increase would act as a sink for and source for ; similar arguments hold in the case of thinning. Divergence affected ice of all thicknesses, and there was a source term for open water. He assumed that depended on the random short-term strain in the ice. By using different values of d and , he was able to show the effects of different processes on g(h). The following is a brief summary of his findings:This study considerably improved our understanding of , but left the following key issues open:The theoretical investigation of the evolution of g(h) was complemented by observations of thickness in the central Arctic, which revealed that for thick ice. Thorndike [16] thus constructed simpler models for the thermal and mechanical processes to explain the observed exponential tail. For the thermal process, he assumed , where is the time scale required to reach . The rate of formation of open water and ridges was assumed to be r. Using dimensional arguments he related H to by:where is some function of F / r. Thorndike [16] argued that because there are a large number of interacting floes, the larger the fraction of a certain thickness the larger the probability of participation in ridging to produce thicker ice. From this logic he arrived at the following form for ;where is the source of open water, is the sink term for the ice that is used for ridging, and the convolution term represents the sum of all interactions that produce ice of thickness h. While this approach overcame limitation 2 from the previous study, the nonlinear integro-differential equation could only be solved numerically. The solutions displayed exponential tails, showing that the simple rules for the thermal and mechanical interactions were sufficient to obtain g(h) for thick ice that were in qualitative agreement with the observations.

When , . Here is the Dirac-delta function and is the “equilibrium” thickness. For a typical profile, g(h) attains the maximum value at .

Choosing and leads to a spread in g(h) on both sides of the maximum, but for and the spread is only on the thinner side.

In order to obtain a steady solution, it is necessary for when . Thus, the solution in this case has very little thin ice.

The solutions with and qualitatively resemble the observed g(h).

A closed form of was still lacking, which prohibited any systematic mathematical analysis of Eq. 2.

It was assumed that ice only from a particular range of thickness could ridge to produce thicker ice, but this is generally not the case [19].

It was difficult to use this framework to study seasonal changes in g(h).

Recently, Godlovitch et al. [5] generalized Thorndike’s approach (Eq. 5) using Smoluchowski coagulation models. These models describe the evolution of a population of particles that can interact in pairs to change their mass, with the rate of coalescence that depends on their mass. Using this formalism, was represented aswhere is the rate kernel and C(K, t) is introduced to ensure that g(h) is normalized. Numerical solutions to Eq. 6 displayed exponential and quasi-exponential tails for a variety of , indicating that the nature of the ridging process is not sensitive to the choice of rate kernel. However, the choice of is important for quantitative comparison with observations [5].

Finally, the World Climate Research Programme Coupled Model Intercomparison Project Phase 5 (CMIP5) models, which use momentum equations for the ice pack and hence a wide range of constitutive models for the ice rheology, are known to poorly represent the spatial patterns of ice thickness [13]. This highlights the utility of taking a probabilistic approach to understanding the large scale behavior of the ice pack, which is what motivated the original theory of Thorndike et al. [14].

A Statistical Mechanics Based Theory [17]

When studying Eq. 2, it is important to realize that there is a separation of length and time scales over which the mechanical processes (e.g., ridging and rafting) act relative to the evolution of g(h). Observations indicate that a region with a length scale of 100 km or more is required to define g(h) [14], whereas the features that result from ridging and rafting extend over up to a few tens of meters in general [19]. Hence, one could construct a theory that neglects the details of the collisions, but takes their net effect into account to study the geophysical-scale evolution of g(h). This line of reasoning led us to use an analogy with Brownian motion to interpret [17]. Now we describe this approach.

The classical problem of Brownian motion concerns the motion of a pollen grain in water [1, 12]. The collisions with the water molecules effect the motion of the pollen grain. Given the length-scale separation between the pollen grain and the solvent molecules, there are an enormous number of solvent-grain collisions over the time scale of the evolution of the pollen grain. Hence, one does not take the individual collisions into account when describing its motion, but only their (appropriately averaged) net effect.

We view the short length and time scales of individual mechanical processes (ridging and rafting) relative to the overall evolution of g(h, t)1 in direct analogy to the collisions of water molecules with a Brownian particle, and thus write asThus, we interpret the mechanical redistribution of ice thickness as the differential form of the Chapman–Kolmogorov equation, or a Master equation. The transition probabilities per unit time and represent deformation processes changing ice from thickness to h and from h to respectively, and . We Taylor expand Eq. 7 thereby transforming Eq. 2 towhereEquation 8 is a Fokker–Planck-like equation that describes the evolution of the probability density g(h, t). Here, and represent the first and second moments of thickness transition events, which, because of our core framework that the events that change the thickness occur very rapidly relative to the overall changes in g(h, t), are constants.

We nondimensionalize this equation by choosing as the vertical length scale; L as the horizontal length scale; as the velocity scale for the horizontal ice velocity; as the time scale for advection of ice floes; , where is the thermal diffusivity of ice, as the diffusion time scale; and , where is the collisional strain rate, as the relaxation time scale. Hence, the remaining terms have the following scalings: , , and . Maintaining the pre-scaled notation and noting that , Eq. 8 is:where . When is solenoidal in the domain R, Eq. 10 becomeswhere . Equation 11 is a Fokker–Planck equation for g(h, t). In this paper, we discuss the analytical and numerical solutions to Eq. 11 with a particular focus on the climatological evolution of the thickness distribution.

Analytical Solutions

Steady Solution

A unique steady solution to Eq. 11 was obtained in [17] as follows. The thermal growth rate was taken to be the solution to the ideal Stefan problem (see e.g., [20]) viz., , where S is Stefan number defined as where , and are the latent heat of fusion of ice, specific heat of ice at constant pressure and the temperature difference across the ice layer, respectively. Using the boundary conditions the steady state solution iswhere and . The prefactor, , is the normalization constant with the Euler gamma function. Finally, we note that for given positive values of q and H, this solution is unique.

Setting to zero the first derivative of Eq. 12 with respect to h yields asthus also providing a derivation of Thorndike’s dimensionless function (Eq. 4) from our steady state solution.

Time-Dependent Solution

As a first step in understanding how the thickness distribution is driven by climatological forcing, we introduce a simple model for the growth rate, with and during growth and melt seasons respectively. We emphasize that this model is intended to be pedagogical, as it does not accurately model the winter growth rate, which depends on the thickness. We use Chandrasekhar’s method [1] to first obtain the fundamental solution to Eq. 11. The method involves computing the characteristics of the advective part of the equation, in which is now a constant, along which we write the resulting–diffusion–equation (see Appendix 1);where and , and the boundary conditions are . The Green’s function that satisfies these conditions isand thus, the time-dependent solution for a given initial condition is given byFor , the solution in terms of h and t iswhere , , , , , and the error function (complimentary error function) is (). Importantly, Eq. 17 demonstrates that even in the case of growth rate being independent of thickness, multiple time scales are generated due to the interaction between thermal and mechanical processes.

Numerical Solutions

Seasonality and climate forcing are introduced by coupling Eq. 11 to the one-dimensional thermodynamic model of Eisenman and Wettlaufer (EW09) [3]. We discretize Eq. 11 using the standard second-order finite-difference formulae and it is integrated in time using the semi-implicit Crank–Nicolson scheme. The equation for is:where is the density of ice, is the incoming shortwave radiative flux, , is the incoming longwave radiative flux, , are the turbulent specific and latent heat fluxes at the upper surface, is the controlled flux perturbation at the upper surface (representing greenhouse gas forcing), and is the oceanic heat flux at the bottom surface. A linearized form of the Stefan-Boltzmann law is used for the outgoing longwave radiative flux and is given by , where T(h, t) is the temperature of the upper surface. Ice export is per year and is represented by . Whilst we have neglected the advection term in Eq. (11), we have incorporated the mean effect of advection on ice export in this manner, but it is not the same as incorporating the full ice velocity field. However, ignoring export leads to relatively minor quantitative changes and no qualitative changes to the results presented here. Finally, we note that because S is large for ice, the energy balance across h is global and hence, for example,  Wm, and ,  Wm are equivalent.

Sea Ice Growth Rate

The typical growth rates from the EW09 model for winter and summer are shown in Fig. 1. Clearly, as this is a Stefan problem the growth rate decreases with increasing thickness, due to the fact that growth rate depends on the amount of heat conducted through the ice layer, which decreases with increasing thickness. The growth rates shown here are similar to those obtained from MU71 [10, 14]. The melt rate is constant for all thicknesses except for , which can be attributed to the accelerated melting of thin ice because of the ice-albedo feedback. The feedback is captured here through the h dependence of the albedo, where the characteristic length scale is the inverse of the spectrally averaged Beer’s extinction coefficient  m [10]. This becomes particularly important as the ice cover thins and because our thermodynamic model does not account for the fraction of energy penetrating into the water column and effectively increasing .

Fig. 1

Typical growth rate profiles for winter (solid line) and summer (dashed line)

Evolution of the Mean Thickness

The mean thickness is defined as:in the same manner as the mean values of all quantities X. Figure 2a shows the seasonal behavior of the dimensional . For all values of the seasonal cycle of is qualitatively the same, with the maximum, , at the end of the growth season in early April, and the minimum, , at the end of the melt season in August. This behavior is in general agreement with observations and with solely thermodynamic models [3, 10]. Importantly, however, for greenhouse gas forcing roughly twice that at which the thermodynamic only component of this model transitions from the seasonal ice state to the ice free state (see Fig. 3 of [3]), here we are still in the perennial state. This exhibits the crucial role that ice mechanics plays in maintaining the ice cover by redistributing thin ice to thick ice-far more rapidly than can thermal growth alone. Indeed, it is not until reaches approximately six times the thermodynamic only transition to the ice free state that the annual mean (see Fig. 8).

Fig. 2

Seasonal evolution of the mean thickness and g(h) versus , and  m throughout. Dash-dotted line  Wm,  m,  m; dotted line  Wm,  m,  m; and solid line  Wm,  m,  m

Fig. 8

Effect of increasing on . Circles simulation; solid line . The mean thickness becomes 0.66 m when  Wm

Seasonal and Climatological Changes in g(h)

For the purpose of discussing the changes in g(h), we define ‘thin’ ice as ice of thickness and ‘thick’ ice as ice of thickness . The fraction of thin ice isWhereas increases by about a factor of two from the winter minimum to the summer maximum (Fig. 3), roughly independent of , the minimum and the maximum also increase by roughly a factor of two as increases.

Fig. 3

Seasonal evolution of the thin-ice fraction with increasing greenhouse gas forcing. Dash-dotted line  Wm; dotted line  Wm; and solid line  Wm

Figure 2b–e show g(h) at the middle and end of the growth and melt seasons respectively. For all , as winter progresses and decreases, both thermal growth and mechanical redistribution drive the spread and rightward motion of g(h) to create more thick ice. As increases during the melt season, g(h) contracts towards thinner ice and the skewness increases, both of which are enhanced substantially as increases (Fig. 4).

Fig. 4

Seasonal evolution of g(h) with increasing greenhouse gas forcing. Dash-dotted line  Wm; dotted line  Wm; and solid line  Wm

We make a qualitative comparison of g(h) with the recent satellite observations from CryoSat-2 [7] in Fig. 5. The data has been averaged over the periods shown in Fig. 5a, whereas data from the model in Fig. 5b are for the particular days shown. The observed spreading of g(h) during winter is explained as above within the framework of the theory; ice growth makes thicker the ice that is formed and subsequently deformed, shifting the peak to the right and broadening the distribution, a behavior that is suppressed as increases. It is the continual deformation of thinner ice to make thicker ice that maintains a thicker ice pack than would be predicted by thermodynamic only models.

Fig. 5

Qualitative comparison with CryoSat-2 observations [7]. a g(h) for first-year ice for the year 2010–2011; the data were obtained over the periods indicated. b g(h) for the whole thickness range for  Wm from our model on the particular days as indicated

Albedo

Importantly, once g(h) is known, all thickness dependent moments can be calculated. A quantity of keen interest is the albedo, whose summer mean values are difficult to model because of the concurrent presence of a wide range of ice thicknesses in the basin. We plot the seasonal evolution of the mean albedo as a function of in Fig. 6. For example, when Wm we see that reaches a maximum (0.671) at the end of the growth season, and a minimum (0.652) at the end of the melt season; a seasonal difference in the extreme values of only , but this translates into a large variation in surface heat balance [4]. Figures 3 and 6 show that and are anticorrelated; and a close observation of the plots reveals a phase difference between them, with leading. Importantly, as increases so too does the amplitude of the seasonal cycle, the peak to peak variation of which has a substantial impact on the radiative forcing and hence the ice thickness.

Fig. 6

Seasonal evolution of the mean albedo with increasing greenhouse gas forcing. Dash-dotted line  Wm; dotted line  Wm; and solid line  Wm

Effects of the Surface Radiative Flux Forcing

The effect of on g(h) can be understood by considering Eq. 18. An increase in results in a smaller growth rate during winter and a higher melt rate during summer, leading to an increase in for all seasons (Fig. 3). Figure 7 shows that the mean growth rate shifts downward with increasing , but the curves are not phase shifted. This shift in is associated with the increase in for all seasons. Figure 10 shows the dependence of the seasonal cycle of the mean ice surface temperature . It is seen that as increases, so too does and thus the winter growth rate decreases. Moreover, the time period during which the upper surface ablates increases with . This combination of effects leads to a decrease in the seasonally averaged mean thickness. Figure 8 shows that decreases exponentially with increasing over the range simulated, and thus should vanish monotonically in the absence of some other feedback, in qualitative, but as noted above not quantitative, agreement with solely thermodynamic models [3, 10].

Fig. 7

Seasonal changes in the mean growth rate of ice with increasing greenhouse gas forcing. Dash-dotted line  Wm; dotted line  Wm; and solid line  Wm

Fig. 10

Variation in the mean temperature with seasons for different . Dash-dotted line  Wm; dotted line  Wm; and solid line  Wm

Finally, given the important shift in the transitions of the ice states (perennial to seasonal to ice free) from solely thermodynamic models to the full mechanical and thermodynamic treatment of this theory, the response of the ice pack to a radiative flux perturbation is clearly different between these approaches. We have quantified the relaxation time scales for different initial conditions (see Appendix 3) and find that there is a range of thin ice fractions, –0.6, for which the relaxation time scale of the ice pack is approximately 50% that of thermodynamic only models; viz., 4 years rather than 10 years. For distributions with much more thick ice the response time scales are controlled by mechanical deformation of thick ice and hence can be much longer.

Conclusion

Using concepts and methods from statistical physics we have transformed the theory of the sea ice thickness distribution, g(h), of Thorndike et al. [14] into a solvable Fokker–Planck-like equation. We have solved the new equation both analytically and numerically using different models for the thermodynamic growth rate f to understand the climatological evolution of g(h). In the simplest case, for the growth and melt seasons, we find an analytical solution (Eq. 17). The solution shows that the interaction of thermal and mechanical processes during the evolution of g(h) leads to the generation of multiple time scales, which in turn affect the evolution. Thus, as previously suggested by Thorndike [15], we do in fact find that g(h) and its moments relax on different time scales, which clearly has important geophysical consequences.

A climatological suite of calculations was performed by coupling the Fokker–Planck equation to the thermodynamic model of Eisenman and Wettlaufer [3]. The transient and time averaged g(h) from our model are in good agreement with the recent satellite measurements over the Arctic basin [7, 17]. As in solely thermodynamic models [3, 10], we find that the stationary state has a mean thickness, , reaching a maximum in early April, which is the end of the growth season, and a minimum in early August, which is the end of the melt season. Due to the combined effects of thermodynamics and mechanics, g(h) spreads during the growth season and contracts during the melt season. As greenhouse gas forcing, , increases, this contraction is enhanced, with a larger skewness and a sharper peak at lower thicknesses. However, this model remains in the perennial ice state for approximately twice that at which its thermodynamic component transitions from the seasonal ice state to the ice free state. This exhibits the crucial role that ice mechanics plays in maintaining the ice cover by redistributing thin ice to thick ice; intuitively, doubling the thickness of thin ice by ridging occurs instantaneously [19] relative to doubling it by thermal growth. Clearly, by such a stage the ice-covered fraction of the Arctic Ocean may be vastly smaller than at present. Nonetheless, these physical processes will persist until other effects, such as changing boundary conditions at lower latitudes, take over. For example, although it is not until reaches approximately six times the thermodynamic only transition to the ice free state that the exponential decay of the annual mean, , is reduced to the decay scale of shortwave radiation  m (Fig. 8).

The seasonal behavior of the thin-ice fraction is anticorrelated with the behavior of , which is correlated with the evolution of the mean albedo . The surface radiative flux perturbation impacts g(h) by decreasing the mean growth rate and the seasonally averaged mean thickness , thereby leading to an increase in . Depending on the initial , the relaxation times for g(h) to reach a stationary state starting from an arbitrary initial condition range from 4 to 10 years. Importantly, for a range of thin ice fractions, –0.5, there is a minimum in the relaxation time of 4 years, which is approximately 50% that of thermodynamic only models. For distributions with much more thick ice, the response time scales are controlled by mechanical deformation of thick ice and thus become much longer.

The results presented here demonstrate veracity of using the methods and concepts of statistical mechanics to study the geophysical-scale evolution of Arctic sea ice. As described in the introduction, the CMIP5 models poorly represent the spatial patterns of ice thickness [13]. The concept of the original theory of g(h) due to Thorndike et al. [14] was to avoid the complexities of unknown ice rheologies in the equations of motion for the ice cover, and to produce a climatologically relevant and easily implementable probability density function of this core geophysical scale state variable. However, implementation was difficult because of the intransigence of the redistribution function . Having solved this problem in our theory, we find solutions that are in good agreement with satellite observations. Therefore, using the present treatment for g(h) in climate models should lead to a more realistic representation of Arctic sea ice within them. The thermodynamic component used here [3] reproduces the seasonal cycle of Maykut and Untersteiner [10], which is the starting point for all subsequent simplifications of the thermodynamics used in climate models. Therefore, the implementation of our approach, which captures both the mechanics and the thermodynamics, in comprehensive models should be of interest.