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, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F_0$$\end{document}ΔF0, increases. The mean ice thickness decays exponentially with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F_0$$\end{document}ΔF0, 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 integral gives the fraction of that area (A/R) that contains ice of thicknesses between h 1 and h 2 , 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 u 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 f 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 d = ∇ · u, the solutions displayed the following features: (a) ψ provided a source of open water; (b) ice of thickness less than a certain value h * 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 h * .
For his second calculation, Thorndike formulated the original equation as a Markov process; and by assuming the forms of f and ψ he solved for the steady state. In constructing the matrices of f and ψ he used the following principle: If ice of initial thickness h i grew either by thermal growth or by ridging to a final thickness h f , then the process that led to this increase would act as a sink for g(h = h i ) and source for g(h = h f ); 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 e in the ice. By using different values of d and e, he was able to show the effects of different processes on g(h). The following is a brief summary of his findings: 1. When d = 0 and e = 0, g(h) = δ(h − H eq ). Here δ(x) is the Dirac-delta function and H eq is the "equilibrium" thickness. For a typical profile, g(h) attains the maximum value at h = H eq . 2. Choosing e > 0 and d = 0 leads to a spread in g(h) on both sides of the maximum, but for d > 0 and e = 0 the spread is only on the thinner side. 3. In order to obtain a steady solution, it is necessary for e = 0 when d < 0. Thus, the solution in this case has very little thin ice. 4. The solutions with d = 0 and e > 0 qualitatively resemble the observed g(h).
This study considerably improved our understanding of ψ, but left the following key issues open: 1. A closed form of ψ was still lacking, which prohibited any systematic mathematical analysis of Eq. 2. 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]. 3. It was difficult to use this framework to study seasonal changes in g(h).
The theoretical investigation of the evolution of g(h) was complemented by observations of thickness in the central Arctic, which revealed that g(h) ∼ e −h/H 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 f (h) = F × (H eq − h), where F −1 is the time scale required to reach H eq . The rate of formation of open water and ridges was assumed to be r . Using dimensional arguments he related H to H eq by: where G 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 δ(h) is the source of open water, −2g(h) 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.
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 as where K (h, h ) 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 K (h, h ), indicating that the nature of the ridging process is not sensitive to the choice of rate kernel. However, the choice of K (h, h ) 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 ψ as Thus, we interpret the mechanical redistribution of ice thickness as the differential form of the where Equation 8 is a Fokker-Planck-like equation that describes the evolution of the probability density g(h, t). Here, k 1 and k 2 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 H eq as the vertical length scale; L as the horizontal length scale; U 0 as the velocity scale for the horizontal ice velocity; t m = L/U 0 as the time scale for advection of ice floes; t D = H 2 eq /κ, where κ is the thermal diffusivity of ice, as the diffusion time scale; and t R ∼ 1/γ , whereγ is the collisional strain rate, as the relaxation time scale. Hence, the remaining terms have the following scalings: Maintaining the pre-scaled notation and noting that t R ∼ t m , Eq. 8 is: 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.

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., f = 1/Sh, where S is Stefan number defined as S ≡ L i /c p T where L i , c p and T 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 g(0) = g(∞) = 0 the steady state solution is where q = τ/k 2 S = /k 2 and H = k 2 /k 1 . The prefactor, N (q) = H 1+q (1 + q) −1 , is the normalization constant with (x) 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 h = H eq as thus 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 f = 1 and −1 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 y = h + φ t and s = t, and the boundary conditions are g(y = φ s, s) = g(y = ∞, s) = 0. The Green's function that satisfies these conditions is and thus, the time-dependent solution for a given initial condition g 0 (y 0 ) is given by For g 0 (y 0 ) = N y 0 e −by 0 , the solution in terms of h and t is and the error function (complimentary error function) is erf(X ) (erfc(X )). 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 where ρ i is the density of ice, are the turbulent specific and latent heat fluxes at the upper surface, F 0 is the controlled flux perturbation at the upper surface (representing greenhouse gas forcing), and F B 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 is the temperature of the upper surface. Ice export is 10% per year and is represented by ν 0 h. 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,

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 h < 1, 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 λ = 0.67 m [10]. This becomes particularly important as the ice cover thins and h ≈ λ because our thermodynamic model does not account for the fraction of energy penetrating into the water column and effectively increasing F B .

Evolution of the Mean Thickness
The mean thickness is defined as: in the same manner as the mean values X of all quantities X . Figure 2a shows the seasonal behavior of the dimensional h . For all values of F 0 the seasonal cycle of h is qualitatively the same, with the maximum, h max , at the end of the growth season in early April, and the minimum, h min , 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  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 F 0 reaches approximately six times the thermodynamic only transition to the ice free state that the annual mean h ≈ λ (see Fig. 8).

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 h ≤ 1 and 'thick' ice as ice of thickness h > 1. The fraction of thin ice is Whereas Φ increases by about a factor of two from the winter minimum to the summer maximum (Fig. 3), roughly independent of F 0 , the minimum and the maximum also increase by roughly a factor of two as F 0 increases. Figure 2b-e show g(h) at the middle and end of the growth and melt seasons respectively. For all F 0 , 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 F 0 increases (Fig. 4).
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 F 0 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.

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 F 0 in Fig. 6. For example, when F 0 = 2 Wm −2 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 2.9%, 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 F 0 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.

Effects of the Surface Radiative Flux Forcing
The effect of F 0 on g(h) can be understood by considering Eq. 18. An increase in F 0 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 F 0 , but the curves are not phase shifted. This shift in f is associated with the increase in Φ for all seasons. Figure 10 shows the F 0 dependence of the seasonal cycle of the mean ice surface temperature T . It is seen that as F 0 increases, so too does T and thus the winter growth rate decreases. Moreover, the time period during which the upper surface ablates increases with F 0 . This combination of effects leads to a decrease in h , the seasonally averaged mean thickness. Figure 8 shows that h decreases exponentially with increasing F 0 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]. 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.3-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, f = ±1 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, h , 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, F 0 , 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 F 0 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 F 0 reaches approximately six times the thermodynamic only transition to the ice free state that the exponential decay of the annual mean, h , is reduced to the decay scale of shortwave radiation λ ≈ 0.66 m (Fig. 8).
The seasonal behavior of the thin-ice fraction Φ is anticorrelated with the behavior of h , which is correlated with the evolution of the mean albedo α . The surface radiative flux perturbation F 0 impacts g(h) by decreasing the mean growth rate and the seasonally averaged mean thickness h , 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.3-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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Exact Solution of the Fokker-Planck Equation
The Fokker-Planck equation for the thickness distribution, g(h), for constant thermal growth rate is where φ = k 1 − τ f , with the boundary conditions g(h = 0, t) = g(h = ∞, t) = 0. The characteristics of the advective part of Eq. 21 are Hence, using the transformation y = h + φt and s = t in equation 21 gives with the transformed boundary conditions; g(y = φs, s) = g(y = ∞, s) = 0. Equation 23 is the diffusion equation for g along the characteristics.
We solve Eq. 23 by first finding its Green's function, g(y, s; y 0 ), which satisfies the following expression where δ(X ) is the Dirac-delta function. Following Duffy [2], one can seek the Green's function in the form g(y, s; y 0 ) = G(y, s; y 0 ) + u(y, s), where G(y, s; y 0 ) is the free-space Green's function given by and u(y, s) is a homogeneous function that satisfies the following boundary condition at y = φs; which ensures the boundary condition g(y = φs, s) = 0 is enforced at all s. The function u(y, s) satisfies Eq. 23 and can be shown to be [2] u(y, Once g(y, s; y 0 ) is known, the solution corresponding to any initial condition g 0 (y 0 ) can be calculated as g(y, s) = ∞ 0 g(y, s; y 0 ) g 0 (y 0 ) dy 0 .
To solve Eq. 31 numerically, we discretize it using the standard second-order finite difference formulae [9], and integrate it in time using the semi-implicit Crank-Nicolson method. This particular method is chosen for its stability and accuracy [11]. The finite difference form of Eq. 31 is where h and t are the ice thickness and time steps respectively. Here, g n i corresponds to g(h i = i h, t n = n t) where h i (i = 1, 2, 3, . . .) and t n (n = 1, 2, 3, . . .) are the discrete values of h and t, and similarly for the remaining terms. We rearrange Eq. 32 as where c 1 = t/4 h and c 2 = k 2 t/2 ( h) 2 . Equation 33 represents a tridiagonal system, which can be solved efficiently [9].
To validate the code, we solve Eq. 31 with the growth rate from the ideal Stefan problem [20]. The solution in this case is [17] Thus, the test used is that starting from different initial conditions, this unique steady state solution should be reached. We initialize g(h) using: where k i is the thermal conductivity of ice, and R(x) is the ramp function defined as The dependence of albedo on thickness is modelled using: where α i and α w are the values of albedo for thickest ice and open water, respectively and λ is inverse of the spectrally averaged extinction coefficient for Beer's law [3]. The radiation climatology used to determine the values of Sensitivity of Surface Temperature to F 0 Figure 10 shows the F 0 dependence of the seasonal cycle of the mean ice surface temperature T . It is seen that as F 0 increases, so too does T and thus the winter growth rate decreases. Moreover, the time period during which the upper surface ablates increases with F 0 . This combination of effects leads to a decrease in the annually averaged mean thickness h .

System Relaxation Time Scales
We define the system relaxation time, Λ R , as the time taken for g(h, t) to evolve to a stationary state starting from an arbitrary initial condition. A knowledge of Λ R is important in answering the following question: Given an initial state of the ice pack, if there is a flux perturbation associated with a change in the environment, how quickly does the system forget its initial condition and reach a new stationary state? The variation of Λ R as a function of the initial condition may also help us understand the interaction between thermodynamics and mechanics that drives the system to the new stationary state. One possible way to answer this question would be to compute the relaxation times of the moments, but from Thorndike's [15] and our calculations (Eq. 17 of the main document) it is clear that even for constant thermodynamic growth rates, g(h) and its moments relax on different time scales.
We calculate Λ R as follows. We start with different initial conditions g 0 (h) = N (a)h a e −h/b varying a and b to obtain different values of Φ. We solve which is Eq. (11) of the main document, with f given by  Figure 11 shows Λ R as a function of the thick-ice fraction ( Φ = 1 − Φ) of g 0 (h), and the dashed vertical line is Φ for the final time-averaged g(h), denoted by Φ f . The two curves, representing different values of a, display similar behavior. This shows that Λ R is a function of Φ to leading order. When Φ < Φ f , so that the thick-ice fraction of the initial condition is less than that of the stationary state, Λ R varies between 4 and 8 years. Therefore, there is a range of of thin ice fractions, Φ ≈ 0.3-0.6 for which the relaxation time scale is approximately 50% that of thermodynamic only models, viz., ∼4 years rather than ∼8 years, where there is an efficient mechanical redistribution of ice thickness, which is faster than solely thermodynamic time scales. However, when Φ > Φ f , Λ R is a non-decreasing function of Φ, and it is increasingly difficult for both thermal and mechanical processes to drive the system to the new state.

Appendix 4: Stability Analysis of the Crank-Nicolson Scheme for the Advection-Diffusion Equation
To study the stability properties of the Crank-Nicolson scheme when applied to an advectiondiffusion equation with constant transport coefficients, we perform a von Neumann analysis. Consider the following advection-diffusion equation where u is some quantity being transported, V is the advection speed, and D is the diffusivity. Using central differences for the spatial derivatives and semi-implicit C-N for time integration, On rearrangement, we find the following u n+1 j − u n j = p 1 u n+1 j+1 − u n+1 j−1 + u n j+1 − u n j−1 + p 2 u n+1 j+1 − 2u n+1 j + u n+1 j−1 + u n j+1 − 2u n j + u n j−1 , where p 1 = V t 4 x and p 2 = D t 2( x) 2 . Assuming wave-like solutions we have u n j = u n e i j k x ; u n+1 j = u n+1 e i j k x ; u n j+1 = u n e i ( j+1) k x , etc.
where i = √ −1 and k is the wavenumber. Using these in Eq. 42 and after some algebra we obtain: where G = u n+1 / u n is the amplification factor, a = 1 − 4 p 2 sin 2 k x 2 , b = 2 p 1 sin (k x), and c = 1 + 4 p 2 sin 2 k x 2 . The squared amplitude of G is For stability, we must have |G| ≤ 1. Assuming a = O(1) and c = O(1), we consider the following cases: This is the amplification factor for the pure diffusion equation for C-N scheme, and implies unconditional stability as |G| ≤ 1 for all cases. 2. When b 1, we have which implies numerical instability.
To ensure numerical stability in our simulations, t and h are chosen such that | p 1 | ≤ 0.4 throughout the year.