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

Abstract

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.

1. 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₂), 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.

2. Materials and Methods 2.1. 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(_Ha+_Hg+1)×_Lx×_Lyboxes where_Hais the number of atmospheric boxes in each vertical column,_Hgis the number of ground boxes in each vertical column and_Lx×_Lyis the number of boxes in each horizontal layer. The position of a box is indexed by its discrete position(x,y,z)∈[1:_Lx]×[1:_Ly]×[−_Hg:_Ha]. The layerz=0corresponds 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_{Fxyz,^{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_Qxy(in W);

• Conduction in the ground between box(x,y,z)and box(x,y,z+1)is noted_Dxyz[T](in W). 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_Rxyz[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,_Rxyz[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_Cxyz(inJ·^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:

_{Cxyz T˙xyz}=_Rxyz[T]+_γxyz,z∈[2:_Ha],

_{Cxy1 T˙xy1}=_Rxy1[T]+_γxy1+_Qxy,

_{Cxy0 T˙xy0}=_Rxy0[T]−_Qxy+_Dxy−1[T],

_{Cxyz T˙xyz}=−_Dxyz[T]+_Dxy(z−1)[T],z∈[−1:−_Hg],

where we have introduced the convergence of atmospheric energy fluxes_γxyz=∑^{x′ y′ z′}(_{F^{x′ y′ z′},xyz}−_{Fxyz,^{x′ y′ z′}}). In usual energy balance models,_{Fxyz,^{x′ y′ z′}}and_Qxy 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_γxyzand_Qxyby maximizing the entropy production created by the turbulent sensible heat exchanges over all the domain and timet∈[0:T]:

σ=∫0T∑x,y,z∈[2:_Ha]_{γxyz Txyz}⏟upperatmosphere+∑x,y_γxy1+_QxyTxy1⏟boundarylayer−∑x,y_{Qxy Txy0}⏟surfacedt,

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:

∑x,y,z_γxyz=0∀t∈[0:T].

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:

A[_γxyz,_Qxy,β]=∫0T∑x,y,z∈[2:_Ha]_{γxyz Txyz}+∑x,y_γxy1+_QxyTxy1−∑x,y_{Qxy Txy0}−β∑x,y,z∈[1:_Ha]_γxyzdt,

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:

A[_Txyz,_T˙xyz,β]=∫0TL[_Txyz,_T˙xyz,β]dt,where

L[_Txyz,_T˙xyz,β]=∑x,y,z∈[1:_Ha]_{Cxyz T˙xyz}−_Rxyz[T]1_Txyz−β+∑x,y_{Cxy0 T˙xy0}−_Rxy0[T]−_Dxy−1[T]1_Txy0−β

is the Lagrangian of the optimization problem. It can be expressed in a more compact form

L[_Txyz,_T˙xyz,β]=∑x,y,z∈[0:_Ha]_{Cxyz T˙xyz}−_Fxyz[T]1_Txyz−βwhere

_Fxyz[T]=_Rxyz[T]+_Dxy(z−1)[T]−_Dxyz[T]

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_Dxyz=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):

∂L∂_Txyz−ddt∂L∂_T˙xyz=0∀x,y,z∈[0:_Ha]

∂L∂β=0.

These equations have to be complemented by the energy budget equations in the ground forz<0. All in all, we then obtain the following set of equations:

_Cxyzβ˙+∑^{x′ y′ z′}∈[0,_Ha]∂_{F^{x′ y′ z′}}[T]∂_Txyzβ+_Fxyz[T]_Txyz2−∑^{x′ y′ z′}∈[0,_Ha]∂_{F^{x′ y′ z′}}[T]∂_Txyz1_{T^{x′ y′ z′}}=0∀x,y,z∈[0:_Ha],

∑x,y,z∈[0:_Ha]_{Cxyz T˙xyz}−_Fxyz[T]=0,

_{Cxyz T˙xyz}=−_Dxyz[T]+_Dxy(z−1)[T],forz∈[−_Hg:−1].

The only inputs of such a model are the discretization of the domain, the radiative model_Rxyz[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 (i.e.,L[_Txyz,_T˙xyz,β]is replacedL[T,_∂μT,β]whereμ=(t,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_Qxyfor 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 modelR[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.

2.2. 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_Ri,i=1,2. They can exchange energy by a “turbulent” heat fluxF=_γ2=−_γ1that will be optimized with MEP. We note their temperature_Tui,i=1,2. The two bottom boxes represent buffers of energy that store/release energy by conduction_Di,i=1,2to the upper boxes (typically the ground). We note their temperature_Tbi,i=1,2._Ri, F and_Diare expressed in W. We consider a local radiative forcing (i.e., that depends only on the box temperature_Ri(_Tui)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_T0ithat physically represents the radiative equilibrium solution of the problem (i.e., the temperature of the upper box_Tuiin the absence of thermal inertia, no convective fluxes_γi=0and no diffusive heat fluxes_Di=0). The linearization around_T0iassumes 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

_Ri=_riT0i−_Tui,

where_riare extensive quantities (in W·^K−1) proportional to the size of the box i. Conduction is represented by

_Di=_kiTbi−_Tui,

where_kiis also an extensive quantity (in W·^K−1). We consider the homogeneous case where_ri=rand_ki=kfori=1,2 . It can be shown that the Equations (14)–(16) reduce to

_Cuβ˙−(r+k)β−r_T0i+k_TbiTui2=0i=1,2,

∑i=12_{Cu T˙ui}−r(_T0i−_Tui)−k(_Tbi−_Tui)=0,

_{Cb T˙bi}=−k(_Tbi−_Tui)i=1,2.

The casek=0 corresponds to the model presented in [25]. In this case, one has

_Tui=r_T0iCuβ˙−rβ≡α(t)_T0i,

such that_Tuiis directly linked to the forcing_T0i, modulated by the functionα(t)≡r/_Cuβ˙−rβ. This solution does not correspond to a “true” local lag due to the local energy storage because_T0iis 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 casek≠0, we have

_Tui=r_T0i+k_TbiCuβ˙−rβ≡α(t)_T0i+kr_Tbi.

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_Tbiis related to_Tui through the diffusion Equation (21). We expect that_Tuiwill be in phase with the forcing temperature_T0iifk≪r. On the contrary,_Tuiwill be in phase with the buffer temperature_Tbiifr≪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 whenk≪rand with various values of_Cu in Figure 3. We observed that there was almost no lag of the upper boxes temperatures, even if we increased_Cu 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.

If the geometry (i.e., the sizes of the boxes) is fixed, the problem depends on the heat capacities_Cuand_Cb, and the conductivity k. In the context of a cyclic forcing, the radiative forcing is characterized by a typical timeτ, a reference temperature〈_T0〉, an amplitudeΔ_T0and 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

_Nb=_Cbkτ,_Nr=_Curτand_Nk=_Cukτ,

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.τ_Nbis the typical lag of the ground (≃ two weeks to two months).τ_Nris the typical lag of the atmospheric response and the radiative forcing (≪1 year)._Nkcontrols the influence of the bottom on the top by diffusion. Then, the evolution equations for our model can be written as

β˙−1_Nr+1_Nkβ−1_NrT0i+1_NkTbiTui2=0i=1,2,

∑i=12_T˙ui−1_Nr(_T0i−_Tui)−1_Nk(_Tbi−_Tui)=0,

_T˙bi=−1_Nb(_Tbi−_Tui)i=1,2.

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_Nb,_Nrand_Nk. 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

q=τF_Cu=1_Nr(_T01−_Tu1)+1_Nk(_Tb1−_Tu1)−_T˙u1,

that can be tracked as a diagnostic variable. It represents the temperature elevation of an upper box due to an energy inputτF.

3. Results

We focused on the case of an antisymmetric forcing:_T01=300+10sin(2πτt)Kand_T02=300−10sin(2πτt)K, and varied the control parameters_Nb,_Nrand_Nkto 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_Tu1and_Tb1. Indeed,_Tu2and_Tb2are the antisymmetric versions of_Tu1and_Tb1with respect to the lineT=300K.

• Influence of_Nb: lag of ground:

We clearly observe in Figure 4a that_Nb 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_Nbincreases, the amplitude of_Tbdecreases and its lag increases. This is a classical behavior for diffusion. When_Nbis 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_Nb.

• Influence of_Nk: bottom to top diffusion:

On Figure 5a, we observe that_Nkinfluences the lag and amplitude of the upper and bottom boxes temperatures (with respect to the radiative forcing). As_Nkincreases, the lag decreases and the temperature amplitudes increase. That is easily understood since a high value for_Nkrepresents 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_Nk (Figure 5b).

• Influence of_Nr:

Results presented in Figure 6a show that_Nr 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_Nrincreases, 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 productionF(1/_Tu1−1/_Tu2)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_T01=290+10sin(2πτt)Kand_T02=310+3sin(2πτt+π4)Kwith_Nb=0.1,_Nr=0.001and_Nk=0.1 . The temperatures and q are shown on Figure 7a,b. The expected lag between the upper boxes and bottom boxes is of order_Nb.

4. 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. 5. 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.

Supplementary Materials

The following are available online at https://www.mdpi.com/1099-4300/22/9/966/s1, Python code.

Author Contributions

Conceptualization, V.L. and D.P.; methodology, V.L. and D.P.; software, V.L.; validation, V.L., D.P. and B.D.; formal analysis, V.L. and D.P.; investigation, V.L. and D.P.; writing-review and editing, V.L., D.P. and B.D.; supervision, D.P. and B.D.; funding acquisition, D.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We thank the anonymous referees for their useful comments that have helped to improve significantly the article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Numerical Algorithm

Equations (25)-(27) are a set of nonlinear differential equations that can be solved numerically using the following procedure:

1. The equations correspond to a dynamical system of unknownx=(_Ta1,_Ta2,_Tg1,_Tg2,β)that we formally write as

x˙=f(x).

2. We discretize the cycle in T discrete time steps_tk=k/T,k=0,...,T-1. The temporal derivative of a function f at time_tkis expressed using the second order approximation(f(_tk+1)-f(_tk-1))/(_tk+1-_tk-1). Formally, the dynamical system is transformed into an equation of unknownX=(x(_t1),...,x(_tT)):

g(X)=0,

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:

(a)

We consider^X(n), a reference point for the nth iteration.

(b)

We linearize the equation around^X(n):

_gi(^X(n))+∂_gi∂_Xj(^X(n))(_Xj-_Xj(n))+O(^(X-X(n))2)=0.

(c)

Keeping only the first order, we obtain the correction for the next iteration:

_Xi(n+1)=_Xi(n)-_{Gij-1 gj}(^X(n))

where_Gij(X)=∂_gi∂_Xj(X)is supposed to be invertible.

(d)

We reiterate the last steps until convergence.

More precisely, we note

X=(_Ta1(_t1),...,_Ta1(_tT),_Ta2(_t1),...,_Ta2(_tT),_Tg1(_t1),...,_Tg1(_tT),_Tg2(_t1),...,_Tg2(_tT),β(_t1),...,β(_tT))

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_tkby the finite difference.

f˙(_tk)≃f(_tk+1)-f(_tk-1)_tk+1-_tk-1=f(_tk+1)-f(_tk-1)2/T.

We therefore obtain the following equations∀k∈[0:T-1]:

-1_NrT0i(_tk)+1_NkTbi(_tk)_Tui2(_tk)-1_Nr+1_Nkβ(_tk)+Tβ(_tk+1)-β(_tk-1)2=0i=1,2,

T_Tbi(_tk+1)+_Tbi(_tk+1)2+1_Nb(_Tbi(_tk)-_Tui(_tk))=0i=1,2,

∑i=12_T˙ui(_tk)-1_Nr(_T0i(_tk)-_Tui(_tk))-1_Nk(_Tbi(_tk)-_Tui(_tk))=0,

where we use the periodic boundary conditionsf(_t0)=f(_tT). This set of equation is clearly of the formg(X)=0. The only nonlinear term that needs to be linearized to apply the Newton's method is

1_NrT0i+1_NkTbiTui2=3_NrT0i+2_{NkTbi(n)Tui(n)2}-21_NrT0i+1_{NkTbi(n)Tui(n)3Tui}+1_{Nk Tui(n)2Tbi}+O(^{(_Tui-_Tui(n))2},...).

It follows that we can write the linearized evolution equations as

21_NrT0i+1_{NkTbi(n)Tui(n)3}(_tk)_Tai(_tk)-1_{Nk Tui(n)2}(_tk)_Tbi(_tk)-T2β(_tk-1)-(1_Nr+1_Nk)β(_tk)+T2β(_tk+1)=3_NrT0i+2_{NkTbi(n)Tui(n)2},i=1,2,

-1_NbTui(_tk)-T2_Tbi(_tk-1)+1_NbTbi(_tk)+T2_Tbi(_tk+1)=0,i=1,2,

T2(_Tu1(tk+1)-_Tu1(_tk-1)+_Tu2(_tk+1)-_Tu2(_tk-1))+(1_Nr+1_Nk)(_Tu1(_tk)+_Tu2(_tk))-1_Nk(_Tb1(_tk)+_Tb2(_tk))=1_Nr(_T01(_tk)+_T02(_tk))

These equations can be put under the matricial form

MX=S,

where M is a constant square matrix of size5Tand S a constant vector that depends on_Nb,_Nr,_Nk, T and^X(n). M can be decomposed into5×5blocs of sizeT×Tas

_M110_M130_M150_M220_M24M25-I_Nb0_M33000-I_Nb0_M440_M51M52-I_Nk-I_Nk0

where I denotes the identity matrix of size T. We have the following expressions for the blocks:

_M11=2diag1_NrT01+1_{NkTb1(n)Tu1(n)3},_M22=diag1_NrT02+1_{NkTb2(n)Tu2(n)3},

_M13=-diag11_NkTu1(n)2,_M24=-diag11_NkTu2(n)2,

_M15=_M25=Z-(1_Nr+1_Nk),T2,

_M33=_M44=Z1_Nb,T2,

_M51=_M52=Z1_Nr+1_Nk,T2,

where we have introduced the square matrix of size T:

Z(x,y)=xy0...0-y-yxy0...0..................0...0-yxyy0...0-yx.

The second member S can be decomposed into 5 sub vectors of size T:

S=(F(_T01,_Tu1(n),_Tb1(n)),F(_T02,_Tu2(n),_Tb2(n)),0,0,1_Nr(_T01+_T02)),

where

F(x,y,z)=3_Nrx+2_Nkz^y2.

This algorithm is implemented in the Python code given in the Supplementary Materials.