A Maximum Entropy Production Hypothesis for Time Varying Climate Problems: Illustration on a Conceptual Model for the Seasonal Cycle

We investigated the applicability of the maximum entropy production hypothesis to time-varying problems, in particular, the seasonal cycle using a conceptual model. Contrarily to existing models, only the advective part of the energy fluxes is optimized, while conductive energy fluxes that store energy in the ground are represented by a diffusive law. We observed that this distinction between energy fluxes allows for a more realistic response of the system. In particular, a lag is naturally observed for the ground temperature. This study therefore shows that not all energy fluxes should be optimized in energy balance models using the maximum entropy production hypothesis, but only the fast convective (turbulent) part.


Introduction
Time varying problems are ubiquitous in climate modeling. Indeed, the amount of radiative energy coming from the Sun, which is the dominant contribution to the energy budget of the Earth, depends on time. The variations of this radiative forcing occur on different time scales (from hours for the diurnal cycle to several thousand years for the variations of orbital parameters). To predict the time evolutions of relevant variables (temperature, humidity, wind, . . .) around the globe, one must consider how the different components of the climate system respond to this forcing, and how these components interact with each other.
The representation of the diurnal and seasonal variations of surface temperature in climate models depends on many processes. The type of surface (ocean or land with more or less vegetation) influences the sensible and latent energy fluxes between the ground and the atmosphere. For the ocean, the radiation that penetrates over a depth of 70 m under the surface allows the ocean to buffer the radiative energy input during daytime and more (days, weeks or months). This is not the case for land because of its opacity. In addition, change of surface type naturally modifies its heat capacity and the kind of energy transfers that can occur in the medium. For example, the heat capacity is more important for water than for ground, and convection usually does not occur in the ground. Overall, change in surface type has an important impact on the mass (water vapor and CO 2 ), momentum (by friction) and energy (sensible and latent heat) exchanges between ground and atmosphere. Obviously, the water and energy budget at the surface, and so the representation of the turbulent sensible and latent heat fluxes, is of particular importance to describe these phenomena.
Most of global climate models use parameterizations to represent sub-grid processes like turbulent fluxes. For example, one of the most popular representations of surface fluxes relies on semi-empirical relations called aerodynamic bulk formulas. They are based on the application of the similitude theory on the atmospheric boundary layer [1]. Such relations aim to predict the relation between the turbulent fluxes between ground and atmosphere as a function of the wind and differences of humidity and temperature. The applicability of these formulas and the estimations of the empirical coefficients have been the subject of numerous studies [2][3][4][5][6]. The issues that are raised by the use of aerodynamic formulas are the same than for the closures used for the representations of turbulent fluxes in the atmosphere and ocean [7,8]. Namely, the applicability of these formulas to various conditions is still an open question. For example, the coefficients appearing in the simplest sensible heat flux parameterizations depend on the stability of the atmosphere. This issue has been addressed by using conditional formulas [9] or by adding more free parameters to fit experimental data over a wider range of conditions [10]. Global climate models are used to perform many experiments such as sensitivity of surface temperature tests [11], but usually fail to represent the last trends of the seasonal cycle [12].
The question of the phase and amplitude of the ground surface temperature during the seasonal (or diurnal) cycle can also be investigated with energy balance models. Such models can be useful to explore how local processes may cause variability in the phase and amplitude of the seasonal cycle of surface temperature [12]. For example, a model based only on a simplified energy budget and a maximum power closure hypothesis to fix the surface energy flux, has been constructed to explain the difference in sensitivity over ocean and land [13]. The maximum power closure hypothesis takes the form of a constrained variational problem that does not involve free parameters. It belongs to a class of thermodynamic optimality hypotheses used in various contexts [14]. Similarly, some climate models have been constructed on the maximum entropy production hypothesis or "principle" (hereafter MEP) [15][16][17][18] (see also the discussion [19] and the reviews [20,21]). All of them are based on the energy budget equation, and MEP appears as a closure hypothesis to determine the energy fluxes.
The statute of the MEP hypothesis is still not clear. It is argued in [22] that MEP is an inference. Other studies suggest that among several steady-states, the most stable is the one of maximum entropy production. This fact is observed for example in numerical simulation of the oceanic circulation [23] or simple bistable chemical system [24]. In spite of the lack of strong physical understanding, MEP models give correct estimations for the meridional atmospheric energy fluxes and temperature profiles in the stationary state [15].
In the context of time-varying problems, it has been observed in [25] that the MEP hypothesis could not reproduce a lag between the forcing and the response in the case of an antisymmetric forcing, which is quite unexpected for a system with thermal inertia. The discussion [26] suggested that there is a link between the maximum entropy production hypothesis and the minimum mixing time of a nonequilibrium system. If true, this statement implies that the MEP hypothesis should not be used to represent slow dynamics, but only the fast turbulent mixing.
Typically, the seasonal cycle is characterized by the amplitude and phase of the periodic surface temperature (relative to the periodic local solar insolation). The difference between the temperature and local insolation phases is called the lag, and the ratio of the amplitudes is the gain [12]. We propose a conceptual model of four boxes where the two upper boxes are forced radiatively, and energy is exchanged with a turbulent energy flux that is optimized with MEP. These two boxes are connected to buffers (bottom boxes) that roughly represent ground and store/release energy by conduction. Interestingly, the lag of the upper box can be modified by a parameter that represents the heating time of the upper box by diffusion. This behavior was not observed in [25]. It suggests that the MEP hypothesis could be used to construct simple climate models for time varying problems if we take care to separate slow processes (energy storage) and fast processes (turbulent mixing), and apply MEP only to the last ones. The outline of the paper is as follows: In Section 2, we present the MEP hypothesis formalism in energy balance models, and we present our conceptual model with four boxes and introduce the control parameters. In Section 3, we study the influence of the control parameters on the temperatures and the energy fluxes represented with MEP. We discuss the results in Section 4.

Time Dependent MEP-Based Climate Models
A climate model is usually a discretization of atmosphere and ground (ocean or land) into boxes that exchange mass, momentum, and energy. Atmospheric models explicitly represent these exchanges, with parameterizations for the representation of the small scales' motions. Energy balance models are only based on the energy budget equation. We represent the climate system by (H a + H g + 1) × L x × L y boxes where H a is the number of atmospheric boxes in each vertical column, H g is the number of ground boxes in each vertical column and L x × L y is the number of boxes in each horizontal layer. The position of a box is indexed by its discrete position (x, y, z) ∈ [1 : The layer z = 0 corresponds to the thin surface of the ground that is radiatively heated/cooled.
To construct an energy balance model, one must specify the energy content as a function of the thermodynamic variables (temperature, humidity, etc.) and properties of the media (heat capacities). We also have to represent the energy fluxes between different boxes ( Figure 1). Those energy fluxes are usually mainly driven by turbulent (sub-grid) motions and are difficult to represent in practice. We note: • Atmospheric (z > 0) and oceanic (z ≤ 0) energy fluxes F xyz,x y z from box (x, y, z) to (x , y , z ) (in W); • Surface energy fluxes due to sensible and latent heat fluxes between the ground and the atmosphere Q xy (in W); • Conduction in the ground between box (x, y, z) and box (x, y, z In oceans, energy transfers can also occur by turbulent mixing and oceanic circulation. We have neglected the horizontal diffusion in the ground due to the large aspect ratio of the Earth's crust. It is also necessary to model the external radiative forcing R xyz [T] (in W) that is a functional of the temperature T and greenhouse gases concentrations. The radiative forcing is also strongly affected by the surface albedo which depends on the kind of surface (land or ocean, type of soil/vegetation, ice cover, . . .). Physically, R xyz [T] represents the net radiative balance of box (x, y, z) (in W). Boxes represent different kind of material, with a more or less important size. Therefore, their properties such as heat capacities C xyz (in J·K −1 ) vary from a box to another.
Considering only the sensible heat term in the expression of the energy of each box, the energy budget can be written as follows ∀x, y: where we have introduced the convergence of atmospheric energy fluxes γ xyz = ∑ x y z (F x y z ,xyz − F xyz,x y z ). In usual energy balance models, F xyz,x y z and Q xy are related to the gradient of temperature and humidity by parameterizations (sometimes only by diffusion, as in references [27][28][29][30][31]). Alternatively, the maximum entropy production hypothesis consists of determining γ xyz and Q xy by maximizing the entropy production created by the turbulent sensible heat exchanges over all the domain and time t ∈ [0 : T ]: taking into account that the energy fluxes F are, by definition, internal to the atmosphere so the sum of its convergence over the domain is zero at any time: Maximum entropy production is here applied to the convergence of turbulent energy fluxes and we directly compute the temperatures. Physically, it represents an efficient heat mixing that is only constrained by the energy budget (Equations (1)-(4)) and the constraint (Equation (6)). The reason why we do not optimize the entropy production created inside the ground layer is that the heat exchange is not ensured by the fast convective processes but by slow conductive processes. For the same reason, only the variables directly involved in the turbulent entropy production are considered as variables of the optimization problem. Following [25], the solution of the problem is a critical point of the following functional: where β is the function of time that is the continuous set of Lagrange multipliers associated with the constraints (6) for all times. We express this functional in terms of the temperatures using the energy budget Equations (1)-(4) as follows: where is the Lagrangian of the optimization problem. It can be expressed in a more compact form simply represents the convergence of energy fluxes that are not optimized with MEP, namely, the sum of the net radiative energy input and the convergence of conductive energy fluxes. The diffusive part of the energy transport in the atmosphere is negligible compared to the convective part so we consider that D xyz = 0 ∀x, y, z ≥ 0. We observe that the sum in the Lagrangian is performed over all boxes where there is turbulent energy fluxes (optimized with MEP). At this point, it is worth mentioning that MEP makes no distinction on the nature of the energy fluxes that are optimized (surface or atmospheric heat fluxes, sensible or latent, diffusive or convective). The solution of the optimization problem satisfies the following Euler-Lagrange equations (see, for example, [32] for an introduction to variational problems in physics): These equations have to be complemented by the energy budget equations in the ground for z < 0. All in all, we then obtain the following set of equations: The only inputs of such a model are the discretization of the domain, the radiative model R xyz [T] and the representation of the non-turbulent (slow) heat fluxes (here the conduction in the ground).
The differences with stationary MEP are that the radiative forcing depends on time and also that the slow variables (here the ground temperatures) can influence the fast variables through (16).
One notes that it is tempting to write the continuous field theory associated with the MEP hypothesis. Is this case, spatial derivatives of the temperature field can appear in the Lagrangian x, y, z)), so the discrete sums over boxes should be replaced by an integration, and Equations (14)- (16) should be replaced by partial differential equations (PDEs). However, contrary to systems usually described by classical field theories (like elasticity or classical electromagnetism), the system is divided here into several sub-systems (atmosphere, surface, ground) which makes it necessary to use finite difference at least to represent the interfaces between sub-systems (the surface energy flux Q xy for example). Thus, expressing the optimization with discrete space coordinates (i.e., dividing the system into boxes) is more pragmatic, even if less elegant. For more, in the end, the climate modeler will need to discretize the equation to solve it numerically, at least for a sufficiently realistic radiative model R[T]. Thus, the derivation used here is more direct when we are interested in the applications to climate modeling, even if it is less usual when one compares it to the field's usual formalism.

Simplified Case
To analyze the physical relevance of MEP hypothesis and get insights into the effect it has on the represented variables, it is instructive to consider the model with only four boxes shown in Figure 2. It represents two columns of two layers. The upper boxes represent the part of the system that is forced radiatively (typically the troposphere and the thin ground surface) by source terms R i , i = 1, 2. They can exchange energy by a "turbulent" heat flux F = γ 2 = −γ 1 that will be optimized with MEP. We note their temperature T ui , i = 1, 2. The two bottom boxes represent buffers of energy that store/release energy by conduction D i , i = 1, 2 to the upper boxes (typically the ground). We note their temperature T bi , i = 1, 2. R i , F and D i are expressed in W. We consider a local radiative forcing (i.e., that depends only on the box temperature R i (T ui ) and external parameters), which is, of course, over-simplistic for a true atmospheric column, but allow for useful simplifications. Furthermore, one assumes that the radiative forcing can be linearized around a temperature T 0i that physically represents the radiative equilibrium solution of the problem (i.e., the temperature of the upper box T ui in the absence of thermal inertia, no convective fluxes γ i = 0 and no diffusive heat fluxes D i = 0). The linearization around T 0i assumes that the temperatures stay close to the radiative equilibrium solution. In this toy model, the atmosphere is reduced to only one layer (that even includes the surface layer). However, for models with several atmospheric layers and a realistic description of the radiative forcing, the radiative equilibrium temperatures are not the same for all the boxes. The radiative forcing then reads where r i are extensive quantities (in W·K −1 ) proportional to the size of the box i. Conduction is represented by where k i is also an extensive quantity (in W·K −1 ). We consider the homogeneous case where r i = r and k i = k for i = 1, 2. It can be shown that the Equations (14)-(16) reduce to The four-box model.
The case k = 0 corresponds to the model presented in [25]. In this case, one has such that T ui is directly linked to the forcing T 0i , modulated by the function α(t) ≡ √ r/ C uβ − rβ. This solution does not correspond to a "true" local lag due to the local energy storage because √ T 0i is naturally in phase with the local forcing, and α(t) does not depend on i so it represents a global response (independent of the position) of the system at time t. In contrast, in the case k = 0, we have The response is still non-local due to the use of MEP and one still have the global modulation α(t). Yet, we naturally have introduced a "true" lag because T bi is related to T ui through the diffusion Equation (21). We expect that T ui will be in phase with the forcing temperature T 0i if k r. On the contrary, T ui will be in phase with the buffer temperature T bi if r k. For real climatic conditions, we obviously observe an intermediate situation.
To illustrate this point, we show the temperatures and the energy flux during the cycle for a typical case when k r and with various values of C u in Figure 3. We observed that there was almost no lag of the upper boxes temperatures, even if we increased C u by several orders of magnitude as in [25]. That is quite unexpected in a real system since the heat capacity naturally controls the thermal inertia of a system. It is explained here by the fact that the energy flux F which maximizes the entropy production adjusts with the change of the heat capacity such that we do not observe a lag. It means that the energy storage (and hence the induced lag) cannot be represented with the MEP hypothesis alone.  . Evolution of temperatures in the first column T u1 (full lines) and T b1 (dashed lines); (a) and energy flux F (b) during the cycle for C b = 1000, τ = 1000, k = 0.001 r = 1.5, various values of C u and antisymmetric forcing temperatures: T 01 = 300 + 10 sin( 2π τ t) K and T 02 = 300 − 10 sin( 2π τ t) K. The forcing temperature T 01 is shown in black. The temperatures in the second column T u2 , T b2 and T 02 are the antisymmetric versions of T u1 , T b1 and T 01 with respect to the line T = 300 K. We see that C u does not influence the temperatures. Here, there is no lag between the forcing and the temperature response.
If the geometry (i.e., the sizes of the boxes) is fixed, the problem depends on the heat capacities C u and C b , and the conductivity k. In the context of a cyclic forcing, the radiative forcing is characterized by a typical time τ, a reference temperature T 0 , an amplitude ∆T 0 and the constant r. These seven physical quantities depend on four physical units (mass, time, length, temperature). Through dimensional analysis, we can characterize the problem through only three dimensionless quantities. We chose the following dimensionless numbers which respectively represent the dimensionless (divided by the cycle's time τ) inertial thermal time scale by conduction in the bottom box, the inertial thermal time scale by radiative forcing in the upper box and the inertial thermal time scale by conduction in the upper box. These parameters control the relative amplitude and lag between temperatures. τN b is the typical lag of the ground ( two weeks to two months). τN r is the typical lag of the atmospheric response and the radiative forcing ( 1 year). N k controls the influence of the bottom on the top by diffusion. Then, the evolution equations for our model can be written asβ where the time has been scaled by τ. The temperatures have not been scaled (i.e., they are still expressed in K) for simplicity. This set of equation has no explicit solution but it can be solved numerically using the Newton's method (Appendix A). We studied the evolution of the temperatures' gains and lags with the variations of the control parameters N b , N r and N k . These quantities depend mainly on the heat capacities (and so the dimensionless parameters introduced above). It is also useful to track the evolution of the energy flux between the upper boxes. For this, we introduce the quantity that can be tracked as a diagnostic variable. It represents the temperature elevation of an upper box due to an energy input τF.

Results
We focused on the case of an antisymmetric forcing: T 01 = 300 + 10 sin( 2π τ t) K and T 02 = 300 − 10 sin( 2π τ t) K, and varied the control parameters N b , N r and N k to study their influences. In this case, the two columns can be viewed as a crude representation of the two hemispheres of the Earth. Due to the symmetry, we only present the temperatures of the first column T u1 and T b1 . Indeed, T u2 and T b2 are the antisymmetric versions of T u1 and T b1 with respect to the line T = 300 K.
Influence of N b : lag of ground: We clearly observe in Figure 4a that N b controls the lag and amplitude of the bottom temperature (with respect to the upper box). This parameter does not significantly affect the upper box temperature (Figure 4a,b) nor q. As N b increases, the amplitude of T b decreases and its lag increases. This is a classical behavior for diffusion. When N b is finite (i.e., k = 0), a local lag is naturally introduced, giving a more satisfactory temperature response than in [25]. q is not influenced by N b .
(a) (b) Figure 4. Evolutions of temperatures in the first column T u1 (full lines) and T b1 (dashed lines); (a) and q (b) during the cycle for varying N b , N r = 0.001, N k = 0.1 and antisymmetric forcing. The forcing temperature T 01 is shown in black. The temperatures in the second column T u2 , T b2 and T 02 are the antisymmetrics of T u1 , T b1 and T 01 with respect to the line T = 300 K. The curves q(t) for different N b cannot be distinguished.
Influence of N k : bottom to top diffusion: On Figure 5a, we observe that N k influences the lag and amplitude of the upper and bottom boxes temperatures (with respect to the radiative forcing). As N k increases, the lag decreases and the temperature amplitudes increase. That is easily understood since a high value for N k represents a large heating time of the upper boxes by diffusion, so they are mainly following the evolution of the radiative forcing. This parameter therefore controls the influence of the bottom ("buffer") on the upper box by conduction. As explained above, it induces a "true" lag of the upper boxes. q decreases with N k (Figure 5b). Influence of N r : Results presented in Figure 6a show that N r controls the lag and amplitude of the upper boxes' temperatures (with respect to the radiative forcing). This parameter has an impact on q (Figure 6b). As N r increases, the amplitude and the lag of q increase. We observed that the lag of q is not the same as the lag of the upper boxes' temperatures (in other words, there is a lag between the upper boxes' temperatures and the energy flux). As a consequence, the local entropy production F(1/T u1 − 1/T u2 ) can be negative during the cycle. It is important to note that it does not necessarily contradict the 2nd law of nonequilibrium thermodynamics which says that the local entropy production due to diffusive heat exchanges is positive, since the definition of entropy production used here encompasses the convective fluxes. We now consider the case of an asymmetric forcing T 01 = 290 + 10 sin( 2π τ t) K and T 02 = 310 + 3 sin( 2π τ t + π 4 ) K with N b = 0.1, N r = 0.001 and N k = 0.1. The temperatures and q are shown on Figure 7a

Discussion
Our model based on the MEP hypothesis has an over-simplistic representation of the Earth system. It was constructed to test the effects of the inclusion of energy storage by conduction in such models, in particular on the lag in the temperature response. The only aspect that MEP keeps about the dynamics is its ability to mix heat. Dynamical constraints have to be combined with the MEP hypothesis (see, for example, [33]) to obtain more realistic results. As discussed in [19], the relevant constraints depend on the problem, and MEP then appears as a useful closure hypothesis that fixes the remaining unknowns, namely, turbulent convective heat fluxes.
When MEP is used, we loose the precise nature and timescales of the different kinds of energy fluxes that are optimized. This means that MEP must obviously be used with caution. This study suggests that we need to have a separation of scales between the global mixing time due to fast turbulent processes, and the scale of interest (which is the period of the cycle here) in order to apply MEP to time-varying problems. This separation of scale is more or less valid for the atmospheric turbulent mixing, but is not true for the energy storage in the ground by conduction that operates on larger timescales. One notes that the separation into fast processes and slow processes recalls the quasi-static approximation, and that this separation depends on the problem considered.
An important question is why entropy production is the relevant quantity to optimize (and why not the power?). This study does not provide a definitive answer. However, one can give the following heuristic explanation: The maximum power hypothesis is interesting for "smart" systems like living organisms (or devices) that are adapting themselves (or adapted) in order to maximize their performances (for example, the energy consumption) despite a reduced efficiency. In such systems, work (the output power) is extracted from the system. On the contrary, convection is not a "smart" adaptative (or adapted) process, but mixes energy (and other quantities) inside the atmosphere so the entropy production is certainly the most pertinent quantity to optimize in this case.
As a final remark, one notes that the maximization of the global entropy production due to turbulent convection implicitly assumes that the turbulent convection is fast enough to mix the atmosphere globally so the global entropy production can be optimized.
We therefore think that it is possible to use MEP to represent turbulent atmospheric motions in more realistic energy balance models (i.e., with a more precise description of the climatic system and its associate radiative budget) to study the seasonal cycle.

Conclusions
We have used the MEP hypothesis in the context of a time-varying problem for a conceptual model. Contrarily to the existing MEP based models, only some energy fluxes were optimized with MEP (the convective ones) and energy was buffered via conduction. This feature naturally led to a lag in the local temperature response, even in the case of an antisymmetric forcing, where no significant lag was observed if we optimized all energy fluxes with MEP. This result suggests that MEP could be used to construct simple climate models for time-varying problems if slow processes, such as conduction in the ground, are decoupled from fast (convective) processes that are fixed by the optimization of entropy production.
We discretize the cycle in T discrete time steps t k = k/T, k = 0, . . . , T − 1. The temporal derivative of a function f at time t k is expressed using the second order approximation ( f (t k+1 ) − f (t k−1 ))/(t k+1 − t k−1 ). Formally, the dynamical system is transformed into an equation of unknown X = (x(t 1 ), . . . , x(t T )): where g is a function that depends on the discretization scheme for the temporal derivative. The problem therefore reduces to the search of zero of the vectorial function g that is solved using Newton's method: We consider X (n) , a reference point for the nth iteration.
We linearize the equation around X (n) : Keeping only the first order, we obtain the correction for the next iteration: We reiterate the last steps until convergence.
More precisely, we note the unknown vector that contains the the concatenated vectors of temperatures and β at discrete times. We approximate the temporal derivative of any function f at discrete time t k by the finite difference.
We therefore obtain the following equations ∀k ∈ [0 : T − 1]: T T bi (t k+1 ) + T bi (t k+1 ) 2 where we use the periodic boundary conditions f (t 0 ) = f (t T ). This set of equation is clearly of the form g(X) = 0. The only nonlinear term that needs to be linearized to apply the Newton's method is 1 N r T 0i + 1 It follows that we can write the linearized evolution equations as T 2 (T u1(t k+1 ) − T u1 (t k−1 ) + T u2 (t k+1 ) − T u2 (t k−1 )) + ( 1 N r + 1 N k )(T u1 (t k ) + T u2 (t k )) − 1 N k (T b1 (t k ) + T b2 (t k )) = 1 N r (T 01 (t k ) + T 02 (t k )) These equations can be put under the matricial form where M is a constant square matrix of size 5T and S a constant vector that depends on N b , N r , N k , T and X (n) . M can be decomposed into 5 × 5 blocs of size T × T as where I denotes the identity matrix of size T. We have the following expressions for the blocks: