Entropy production and multiple equilibria: the case of the ice-albedo feedback

Nonlinear feedbacks in the Earth System provide mechanisms that can prove very useful in understanding complex dynamics with relatively simple concepts. For example, the temperature and the ice cover of the planet are linked in a positive feedback which gives birth to multiple equilibria for some values of the solar constant: fully ice-covered Earth, ice-free Earth and an intermediate unstable solution. In this study, we show an analogy between a classical dynamical system approach to this problem and a Maximum Entropy Production (MEP) principle view, and we suggest a glimpse on how to reconcile MEP with the time evolution of a variable. It enables us in particular to resolve the question of the stability of the entropy production maxima. We also compare the surface heat flux obtained with MEP and with the bulk-aerodynamic formula.


Introduction
A very broad class of problems in climate modelling consists of studying the evolution of a particular field (e.g. surface temperature, precipitation,etc) when an external parameter, or forcing, is varied. Most of the time, the response to this variation is not linear. Feedbacks can amplify or damp the effect of the initial perturbation. One of these feedbacks aroused a proficient branch in scientific literature in the 70s', when Budyko and Sellers simultaneously suggested that the interaction between sea ice and climate could have dramatic consequences. Indeed, the higher the global temperature on Earth, the less the ice cover is likely to extend, and thus the lower the albedo. A lower albedo leads in turn to a higher global temperature, and so on and so forth until all the ice is melted. Stimulated by this pioneer work, the questions of the stability of the climate as well as the consequences such feedbacks might have for understanding paleoclimates were extensively studied, using the whole hierarchy of models, from the most simple Energy Balance Models (EBMs) to the complex General Circulation Models (GCMs).
Using 1D EBMs, Budyko and Sellers had found two stable equilibrium positions for the edge of the ice cover, one corresponding to the present climate and one to a fully ice-covered Earth [3,55]. A large part of the subsequent work was concerned with verifying that these results still held with various different versions of the Budyko-Sellers models, with different heat transport parameterizations, temperature dependance expressions in the planetary albedo, numerical schemes,... [10, 54, 23, 39, 17, e.g.]. Some elegant analytical solutions were found for these models [5,39,40], and various mathematical methods were applied to determine the stability of the equilibria [19,58,16,4,42]. Owing to the extreme sensitivity of climate to variations in the solar constant found by the first studies, the precise position of the tipping point between present climate and a deep freezed Earth was of primary concern. Further investigation by [29] and [43] revealed that the sensitivity was much less than initially thought. A fundamental question raised by these results was that of the transitivity of the climate system in Lorenz's terminology [30,31], and the difference between forced and free fluctuations [54,19,14]. For a comprehensive review of the various models, parameterizations and problems pertaining to Energy Balance Models and the ice-albedo feedback, the reader is referred to [41].
In this contribution, we will first give a brief account of the reformulation of these questions with the vocabulary of dynamical system theory: how do multiple equilibria arise from the ice-albedo feedback, what does the bifurcation diagram look like, etc. The model used here is a two box energy balance model with a simplified radiative transfer using the Net Exchange Formulation (see e.g. [9]), and a bulk aerodynamic formula for the surface heat flux. In a second step, we draw an analogy between this dynamical system view and the results obtained when predicting the surface heat flux from the Maximum Entropy Production (MEP) principle. The MEP principle, as originally expressed by [46,47,48] for the climate system, provides a variational principle to compute energy fluxes that are not otherwise constrained by the laws of physics. Originally, Paltridge and others applied MEP to the meridional energy transport [46, 47, 20, 18, 32, e.g.], but other studies [44,52] indicate that it may be valid on the vertical also.
As noticed by [43], [6] and [14], the bifurcation giving birth to multiple equilibria in the case of the ice-albedo feedback has a fundamentally radiative nature, and has nothing to do with transport properties of the atmosphere. This encourages one in thinking that a zero-dimensionnal model is sufficient to capture the structure of the mechanism while avoiding the use of more cumbersome mathematics (namely the Sturm-Liouville theory, required for one-dimensional models such as [19]). Therefore we will restrict ourselves here to this idealized case. Note also that most of our work could be transposed easily to other feedbacks, like the water-vapour feedback.
2 The ice-albedo feedback, multiple equilibria and the hysteresis cycle: the dynamical system approach

A simple two-layer EBM using the Net Exchange Formulation
We use a slightly different formulation of the model described in [24], as represented in Fig. 1. A grid cell is characterized by a surface temperature T g and an atmospheric temperature T a , and we note Ψ SW gs (respectively Ψ SW as ) the flux of solar energy received by the ground (respectively absorbed by the atmosphere). Radiative exchange use the Net Exchange Formulation, in which the basic objects are not energy fluxes at a given level but rather the energy exchange rate between two layers in the atmosphere or between one layer and a boundary surface (see [9]). Ψ IR ag is the net energy exchange rate between the ground and the atmospheric column per unit surface (ie the greenhouse effect), and Ψ IR sa (respectively Ψ IR sg ) is the cooling to space term for the atmosphere (respectively the surface). The net energy exchange rates per unit surface are expressed as functions of T g and T a as: where σ is the Stefan-Boltzmann constant, α is the surface albedo, t, s, s * ,s are radiative coefficients, S is the solar constant, ξ accounts for the annual mean zenith angle of the sun and µ is the Elsasser factor (see [24] for a derivation of the equations and a discussion of the coefficients).
In addition to radiation, energy is exchanged due to atmospheric and oceanic transport as well as surface heat fluxes. Let us merge all these energy transfer modes into two variables: γ a (respectively γ g ) represents the net convergence (the opposite of the divergence) of energy into the atmospheric cell (respectively the surface layer). Writing ζ a for the atmospheric convergence (this variable was designated by ζ in [24]), ζ o for the oceanic convergence (this was not taken into account in [24]), and q for the surface heat flux, we have Knowing the convergence of energy in each cell -atmosphere or groundit is in general not possible without further assumptions to separate the contribution due to surface fluxes, atmospheric transport, and oceanic transport when applicable. Of course, over land, it is reasonable to assume that γ g is just the surface energy flux (ie ζ o = 0) , and then γ a + γ g is the convergence of energy due to the atmospheric flow. In this study, as we will only use the zero-dimensional version of this model, we will always have ζ a = ζ o = 0, and thus γ a = −γ g = q.
At steady-state, the energy balance equations for the atmosphere and the surface read R a (T a , T g ) + γ a = 0, where and In this form, the steady-state equations (8)-(9) cannot be solved since γ a and γ g are unkown. In the next section we introduce a parameterization of these quantities as functions of T a and T g . In section 3, we use the MEP principle to compute them.

The zero-dimensional model with bulk aerodynamic formula
In the case of a zero-dimensional, two-layer model considered here, the net convergence of energy in the atmospheric box (ie the divergence of the diabatic heating at the surface, γ a = q = −γ g ) can be simply interpreted as the surface heat flux. In this section, we adopt a bulk aerodynamic formula [49] to express this flux as a function of the temperatures T a and T g : where C D is the drag coefficient, u s is the surface wind speed and c pa is the heat capacity of the atmosphere per unit surface area (similarly c pg is the heat capacity of the ground). Now the model can be seen as a two-dimensional dynamical system: with and Our main interest here is to find the equilibrium positions of the system, ie the fixed points of the dynamical system, given by the roots of F , and to study their stability. Of course, the dynamics of a two-dimensional dynamical system can be more complex than just a relaxation to an equilibrium position (although it is still rather gentle, see [22] for example), contrary to onedimensional dynamical systems. Still, let us note here that the first equation is a fixed point of the system. Thus the number of fixed points of the two-dimensional system is exactly the number of roots of the scalar equation For simplicity, we will consider here the projection of the dynamical system (13) onto the T g axis:Ṫ As just explained, this dynamical system, although not mathematically equivalent to the full dynamical system (14), has the same equilibrium positions. Physically, this simplification is motivated by the fact that the atmosphere can be assumed to reach equilibrium very quickly, hence the evolution of T a is enslaved by the dynamics of T g . In other words, the system (17) is just the system (14) with c pa = 0.

Multiple equilibria
The values of the coefficients used here are reproduced in table 1. Taking for the albedo the fixed value α 0 = 0.15, the system only has one fixed point, as plotting the function F 2 (f (T g ), T g ) = 0 clearly shows (see Fig. 2). In this case, the equilibrium is at a global mean surface temperature of T 0 g ≈ 288K. But in reality, the higher the global mean temperature, the lower the extent of the regions that can sustain an ice-cover. This positive feedback can be encoded in the following temperature dependance for the albedo : where α F (respectively α I ) represents the value of the planetary albedo over an ice-free (respectively fully ice-covered) area, and T 0 and ∆T are parameters determining the transition from ice-free to ice-covered conditions (see Fig. 3). One could simply use a step function between ice-free and ice-covered albedo values, or a piecewise linear function, but we choose this expression because it depends smoothly on the temperature.
Replacing α in Eq. (14) with (18) yields a new dynamical system where the fixed points are again determined by the conditions, g being defined similarly to f (or obtained by substitution of the albedo function into f ), Plotting the curve G 2 (g(T g ), T g ) as a function of T g (Fig. 4) shows that for certain values of the solar constant, three solutions coexist. This range can be determined to be approximately 0.98S 0 ≤ S ≤ 1.08S 0 . Outside this range, only one solution subsists. For the present value of the solar constant, S = S 0 , for instance, these equilibria correspond to a fully glaciated Earth (snowball state) T S g ≈ 249K, an ice-free Earth T P g ≈ 287K which can be identified with the present climate, and an intermediate glacial state For a low value of the solar constant (e.g. 0.95S 0 ), only the snowball state T S g subsists. Similarly, at high solar constant (e.g. 1.1S 0 ), the only equilibrium is found on the ice-free branch T P g . A fixed point X * of the dynamical systemẊ = F (X) is said to be (linearly) stable if all the eigenvalues of the jacobian of F are negative (see [1] for a complete classification of the two-dimensional fixed points). In this model we find that T P g and T S g are always stable nodes when they exist, while T G g is a saddle-point. The stability can also be read directly on Fig. 4 for the 1D-reduced system: stable equilibria correspond to roots of the function with negative derivative, while at the unstable equilibrium, the curve crosses the x-axis with an upward slope.
Summarizing the above results, Fig. 5 represents the curve of the fixed points when sweeping a large range of values for S : it is the bifurcation diagram of the dynamical system. Creation of a pair of stable/unstable equilibria at the tipping points 0.98S 0 and 1.08S 0 is called a saddle-node bifurcation. Thus the hysteresis curve obtained for the temperature stems from the bifurcation structure of the dynamical system as two back-to-back saddle-node bifurcations. It is noteworthy that this figure does not depend upon the particular coefficients choice in the bulk formula, nor on the greenhouse effect. Would we set q baf = 0 (radiative equilibrium with greenhouse effect) or/and t = 0 (greenhouse effect shut down), the hysteresis curve would remain qualitatively the same.

Potential for the dynamical system
The full two-variables dynamical system (13) cannot be expressed as the gradient of a potential function, but its one-dimensional projection can, like any other one-dimensional dynamical system. Let us thus introduce the potential V (defined up to an additive constant) such thaṫ Fixed points of this dynamical system correspond to critical points (ie extrema in this 1D case) of the potential. The stability criterion becomes that stable fixed points are minima of the potential: while its maxima are unstable fixed points. Figure 6 shows the shape of the potential for different values of the solar constant. At low solar constant (e.g. 0.95S 0 ), the potential has only one critical point, a minimum at T ≈ 245K. Increasing the value of the solar constant levels down the potential curve, until a second local minimum appears (along with a local maximum) with T above the freezing point, around S ≈ 0.98S 0 . At S = S 0 , it is clear that the potential has two minima at T ≈ 250K and T ≈ 290K and a maximum at T ≈ 275K. Further increase of the solar constant leads to a deeper minimum at T > 0˚C while the minimum at T < 0˚C becomes shallower. Around S ≈ 1.08S 0 , the minimum at T < 0˚C disappears (it annihilates with the local maximum) ; for S = 1.1S 0 , the only minimum is found at T ≈ 300K.
Note that, as expected, the critical points of the potential obtained for the three values of the solar constant considered here match with the values of Fig. 4. Also, the number of critical points of the potential changes at the bifurcation points of the dynamical system.

The entropy production rate and the ice-albedo feedback
In this section, we do not use anymore the bulk aerodynamic formula for the surface flux γ a = −γ g , but the Maximum Entropy Production Principle, as described in [24]. The first application of the MEP principle to climate is found in [46], where the meridional energy transport in a zonally averaged box-model is chosen so as to maximize the entropy production. The resulting climate is in striking accordance with observations. In spite of successful applications in other areas as well, the domain of validity of the MEP principle remains unclear due to the lack of a fully convincing proof (see [7,8] and the comments in [21,2]). More details on theoretical issues and practical use can be found in [45,27,35].

The entropy production rate in zero-dimensions
Let us consider the model of Sect. 2.1 and introduce the entropy production rate per unit surface Substituting Eqs. (8)-(9) into Eq. (24) for γ a and γ g , σ can be considered as a functional of the temperature field. We are looking for its maxima subject to the constraint γ a + γ g = 0.
The sum of equations (8) and (9) can thus be solved for T a as a function of T g , and the entropy production rate σ is simply a function of one variable. Its graphical representation for the set of parameters given in table 1 (fixed albedo α 0 ) is shown in Fig. 7. It is clear that there is only one local maximum, corresponding to a surface temperature T g ≈ 295K.
Now, replacing in the equations the constant albedo α 0 by the temperaturedependent albedo (18), the resulting entropy production rate curve is plotted in Fig. 8 for different values of the solar constant.
Unlike the potential for the dynamical system in Sect. 2.4, the entropy production rate always has at least two local maxima and a local minimum. In fact, over a rather narrow range, estimated to be 0.95S 0 ≤ S ≤ 1.005S 0 , the entropy production rate even has three maxima and two minima. This is even clearer on the contour plot of the entropy production rate as a function of T g and S/S 0 (Fig. 9). Hence, there is indeed an analogue of the fold of the potential in the classical dynamical system picture in the context of the entropy production surface, but the values at which it takes place do not exactly correspond.
Besides, a large portion of the curve on Fig. 8 lies under the abscissa axis: for the corresponding range of temperature values, the entropy production rate is negative, contrary to what the second law of thermodynamics states (or more precisely its extension to non-equilibrium systems). It seems reasonable to impose the condition thereby restricting the range of values T g can actually take. In this case, this is equivalent to requiring that the surface heat flux goes from hot to cold. With this additional constraint, the range of possible values of the solar constant allowing for coexistence of multiple equilibria (two or three) can be determined approximately: 0.8S 0 ≤ S ≤ 1.12S 0 .

Stability of the MEP states
In the classical understanding of the MEP principle, the rate of entropy production σ is a function defined on the manifold of steady-states which reaches a maximum at the most probable state. In the presence of several local maxima, it is generally believed that the final equilibrium state of the system will be the global maximum. In our case, there are three local maxima of the entropy production rate for the present value of the solar constant, as discussed in the previous section. We know from the dynamical system approach that there can indeed be several steady-states (that coincide with the positions of the entropy production maxima, as discussed in the previous section) for a given set of parameters, and the actual steady-state of the system is determined from the initial conditions. In the absence of fluctuations, the system remains in this state. Hence it is certainly not sufficient to retain the global maximum of the entropy production rate as representing the final state of the system. Instead one must find a practical way to select a local maximum for given initial conditions. As a particular case, we would obtain a way to distinguish between local entropy production maxima representing dynamically stable steady-states and dynamically unstable ones.
This involves the introduction of time in the MEP formulation. So far, there was no mention of time in the MEP approach as we were only concerned with steady-states. Even though we claim that the entropy production maxima correspond to equilibrium points, −σ is by no means a potential for the dynamical system. Indeed, the dynamics of the system is simply given by the first law of thermodynamics. Here, it reads Similarly to the steady-state entropy production rate, we can define the instantaneous entropy production rate: using Eqs. (27)- (28). Note that the instantaneous entropy production rate σ i and the steady-state entropy production rate σ coincide at steady-state.
As σ i appears as the natural generalization of σ taking into account the time derivative of the dynamical variables, we suggest that the system may follow the trajectory maximizing the instantaneous entropy production rate, seen as a function of the time-dependent unknown fluxes γ a , γ g (always linked by the relation γ a + γ g = 0). This approach is very similar to what [25] advocates for.
In practice, it is easier to reformulate the above suggestion with a timediscretized system (see Fig. 10). Let us consider two snapshots of the system separated by a finite time interval dt. We note T t a , T t g the values of the air and surface temperature at time t. The instantaneous entropy production rate becomes: Suppose we know the state of the system at time t − 1 (ie T t−1 a and T t−1 g are given). Then our postulate is that T t a and T t g can be chosen so as to maximize σ t i (with fixed T t−1 a and T t−1 g ) subject to the constraint γ t a + γ t g = 0. Iterating this process leads to a trajectory maximizing the instantaneous entropy production rate at each timestep, starting from a given initial condition.
Integrating the system with this method, initialized in the vicinity of the different maxima of the entropy production rate at steady-state, provides a criterion for stability: it is found here that the warm branch as well as the snowball branch of Fig. 9 are stable, while the intermediate branch is unstable. The maxima of the entropy production and their stability are plotted as functions of the solar constant on Fig. 11, analogously to Fig. 5. This result draws the final line in the parallel between the dynamical system approach and the MEP approach. Note that the limits of this analogy are reached at some points: figure 11 cannot be considered as a usual bifurcation diagram. As a consequence, the lines of existence of the maxima need not depend continuously on the parameter, and for certain values of the parameter (for example S ≈ 0.9S 0 ), two stable maxima coexist with no unstable manifold to separate them.
The trajectory maximizing the instantaneous entropy production rate in the way explained above thus yields stability properties for the different steady-states that are consistent with the dynamical system approach. Hence, it seems legitimate to use this hypothesis as a relaxation equation, in a similar fashion as [53]. However, there is no certainty that the system actually follows this maximum instantaneous entropy production trajectory. It would be very valuable to investigate the range of validity of this new application of the MEP principle in future studies, theoretically or on other examples. We can already adduce some material to support our relaxation equations approach. In fact, the only novelty as compared to the common use of MEP in the steady-state context is the inclusion of time derivatives of the dynamical variables in the entropy production rate. But one can simply consider these time derivatives as known fluxes, playing exactly the same role as R a (T a , T g ) or R g (T a , T g ). The only difference is that computing these fluxes requires that we consider a bigger system (here simply the state of the system at times t − 1 and t), even though the number of unknowns in the big system remains the same (T t a and T t g , whereas T t−1 a and T t−1 g are fixed). In this respect, there is no fundamental difference between the time dimension and any geometric dimension, which are customarily included in MEP models.
Alternatively, one could consider the total entropy production rate (ie the integral of the instantaneous entropy production rate over time) as a functional of trajectories and claim that the system follows the trajectory that maximizes this functional subject to the relevant constraints ( [12,13], [11] and [36] have developed this idea in the case of Markov chains by maximizing the information entropy as a function of both the probability of the states and that of the transition rates). As we should show in a forthcoming study, this is particularly suitable for periodic phenomena, such as the seasonal cycle. Regarding the stability of the steady-states, we expect this method to yield the same results as the maximum instantaneous entropy production relaxation used here.

Surface heat flux and Snowball Earth deglaciation
In the case of the first section, the surface heat flux is parameterized as a function of T a and T g . As a consequence of this strong constraint, one could draw a bifurcation diagram for q baf very similar to Fig. 5, with relatively weak surface heat flux q S baf for low solar constants (around 20 W.m −2 ), strong surface heat flux q P baf at high solar constants (around 100W.m −2 ), with an unstable branch q G baf linking the two. On the contrary, the surface heat flux obtained through the MEP procedure q mep is much less constrained by the temperature gradient. Figure 12 shows the surface heat flux as a function of the temperature gradient T g − T a for both cases: q baf and q mep . It is clear that the two differ completely, not only because the temperature gradients in the different climates are very different, but also because the shape of q mep as a function of the temperature gradient is far from linear. Note that in the MEP snowball state, although the temperature gradient is relatively high, the surface flux remains very low. On the warm branch for the MEP state, high values of q mep are obtained for high values of the solar constant. Hence, decreasing the solar constant brings the surface flux down, until the point where only the snowball state survives, with a similar low value of the surface heat flux.
This discrepancy between the two graphs is likely to be significant: it has been suggested that the suppression of the vertical temperature gradient in the snowball state numbers amongst the reasons that make deglaciation of the snowball Earth so difficult [50,51,28]. Indeed, the temperature inversion isolates the surface from all the forms of energy exchange: the greenhouse effect can only warm the surface when the air aloft is colder, latent heat plays a very limited role in this very dry atmosphere, and the sensible heat flux is also restricted by the vertical structure of the atmosphere. [50] points out that a crucial role may be played by the surface fluxes parameterization and the convection parameterization. Here the simplicity of the model does not allow us to discuss the static stability, nor to come up with a clear explanation of the questioning Fig. 12, but it does certainly reinforce the idea that surface heat fluxes parameterization can play critical parts on important paleoclimate problems. In the case of the MEP surface heat flux, our results tend to indicate that it would be possible for the snowball earth to withstand a vertical temperature gradient higher than expected with very little loss in the form of sensible heat, thereby damaging the thermal shield of the surface layer mentioned above.
On a similar note, [34] performed a thorough investigation of the thermodynamic properties of the snowball Earth as compared to warm climates in the model of intermediate complexity PLASIM [15], using the formalism of non-equilibrium thermodynamics applied to the climate system as described in [33]. Computation of the thermodynamics efficiency, irreversibility and material entropy production clearly characterizes distinct thermodynamic regimes for the snowball Earth and ice-free climate. Our remarks about the surface heat flux in snowball conditions add up to their thermodynamic analysis.

Conclusion
The analogy developed in this study leads to some enlightening conclusions. First, about the ice-albedo feedback in itself, it provides a variational principle different from those previously suggested, with a thermodynamic motivation. On the contrary, all the candidates for variational formulations of the problem examined previously were rather ad hoc potentials for the dynamical system. The parallel between potentials properly speaking, which fully describe the dynamics of the system, and the entropy production rate, which only characterize equilibrium states, was pushed one step further with the introduction of a method to integrate a trajectory using the MEP principle. In particular we have shown that this method predicts the correct stability for the MEP predicted equilibria. We also investigated the behaviour of the surface heat flux in the snowball state. The results hint that MEP might prove useful in such extreme situations where the usual parameterizations face important difficulties. However, the highly simplified model considered here does not allow us to conclude against or in favour of the MEP parame-terization, as compared to the bulk-aerodynamic formula.
As far as the MEP conjecture is concerned, our work adds up to the relatively short list of efforts up to now (essentially [56,57] and [26]) to sort out how the principle should be understood in the presence of multiple entropy production maxima. [57] suggested that a dynamical system, in their case the thermohaline circulation, when multiple steady-states are available, should move to the most dissipative one. [37,38] showed strong limitations to this interpretation in full generality. Here, we find that steady-states of a system with unknown turbulent fluxes correspond to local maxima of the entropy production seen as a function of the unknown fluxes. The stability of these maxima does not seem to depend on the numeric value of the entropy production at that point. Instead, we suggest that the question of the dynamic stability can be investigated by a relaxation process maximizing the instantaneous entropy production rate.  Figure 1: A grid cell of the model, adapted from [24]. Ψ ν ij are the energy exchange rates per unit surface due to radiative transfer (see text), q is the surface heat flux and ζ a is the atmospheric energy convergence. Over the oceans, there is also an oceanic energy convergence ζ o .   Figure 10: To discuss the stability of the steady states predicted by MEP, we need to extend the principle to obtain a time-dependent formulation. This is done by maximizing the instantaneous entropy production rate. To compute the time derivative of the temperature, we consider it as a known flux in time seen as a geometric dimension of the space upon which MEP operates (see text). In green, the fluxes that can be computed from the state variables (T t−1 a , T t a , T t−1 g , T t g ). In red, the unknown flux obeying MEP.