A radiative–convective model computing precipitation with the maximum entropy production hypothesis

All climate models use parameterizations and tuning in order to be accurate. The different parameterizations and tuning processes are the primary source of difference between models. Because models are tuned with present observations of Earth, they may not accurately simulate climates of other planets or palaeoclimate. A model with no adjustable parameter that happens to fit today's observations is probably more universal and should be more appropriate to model palaeoclimate. However, to our knowledge, such a model does not exist or is yet to be developed. This paper aims to improve a parameter-free radiative–convective model that computes a realistic temperature vertical profile to compute the water cycle, giving a value on average tropical precipitation. Although it is known that the radiative transfer constrains the order of magnitude of precipitation, no parameter-free model has yet been able to compute precipitation. Our model finds a precipitation value closer to observations than similar radiative–convective models or some general circulation models (GCMs).

1 Introduction

Historically, climate models have evolved from elementary conceptual models to energy balance models (EBMs) and then to radiative–convective models (RCMs) and, after that, general circulation models (GCMs) and, finally, the state-of-the-art Earth system models (ESMs) (Edwards, 2011), with a constant increase in complexity and calculation rate.

Researchers generally consider GCMs and ESMs “the best” models because they account for many phenomena. Indeed, if some specific part of the model does not fit well enough with observations, it is always possible to spend time adding more complexity and making it fit better. There is a hope that if we put enough work into it, GCMs or ESMs will be very close to observations. Furthermore, these models cover the entire Earth, accurately describing the position of oceans and continents, the orography, the cryosphere, and the vegetation with a resolution now below a hundred kilometres. Thus, they are beneficial in answering specific questions, for example, how crop yield would evolve in a particular area within this century. Today's ESMs predict very robust temperature changes for increasing levels of CO₂ (see Lee et al., 2021 Fig. 4.19) for most regions of the globe. On the contrary, they do not predict robust changes in precipitation (Lee et al., 2021 Fig. 4.24). Individual models show opposite signs of precipitation changes in some regions (Fig. 4.42d of Lee et al., 2021). Even when looking for global mean changes when increasing CO₂, temperature shows much less uncertainty than precipitation (Lee et al., 2021 Fig. 4.2a–b). The reason for such uncertainty lies in the difficulty of parameterizing the equations for water fluxes.

Indeed, the atmospheric part of GCMs or ESMs is based on the Navier–Stokes equations, whose length scales range from L=10³ km to less than η=10 mm, the viscous or the Kolmogorov scale of the atmosphere. The number of modes required to model every scale is $N=(L/\mathit{\eta }{)}^{\mathrm{3}}\approx {\mathrm{10}}^{\mathrm{24}}$, by far unreachable by today's computers (we would need 10⁸ times the full storage capacity of today's supercomputer to store one time step). To deal with this problem, climate computers only integrate large scales with the Navier–Stokes equation, and sub-grid processes are parameterized differently for every model, leading to different results, although the physics of every model is the same. How much the parameterizations affect our ability to truly predict climate is a deep and open question (Hourdin et al., 2017; Dommenget and Rezny, 2018). Today's ESMs have hundreds of adjustable parameters.

Looking at global mean evaporation in the tropics, GCMs from the atmospheric model intercomparison project have a mean positive bias of 20 W m⁻² (Zhou et al., 2020). Radiative–convective models compute similar values of mean evaporation to GCMs (Betts and Ridgway, 1988; Rennó et al., 1994; Takahashi, 2009). While radiative–convective models approximate spatial mean quantities well, they are sometimes close to observations at finer scales (Jakob et al., 2019). O’Gorman and Schneider (2008) showed that a GCM gives similar output to a radiative–convective model (using the same simple radiative code and adequate parameterization of surface wind speed for the radiative–convective model); for example, a similar relationship between temperature and precipitation is found (Figs. 3 and 4 of O’Gorman and Schneider, 2008). Equation (19) of Takahashi (2009) gives an upper bound for the surface latent heat flux (and so to evaporation, equal to precipitation in a stationary setting). The exact evaporation value depends on the sensible heat flux and surface radiative transfer. However, if the same drag coefficient and wind speed are used to compute surface sensible and latent heat flux, then because of the exponential dependence of mixing ratio with temperature, latent heat flux dominates sensible heat flux (see Eqs. 7 and 8 of Rennó et al., 1994). Thus, only the knowledge of the surface radiative transfer gives an upper limit to the evaporation. Therefore, in our goal to compute precipitation with a parameter-free model, it is relevant to use a radiative–convective model.

The efficiency of radiative–convective models in reproducing accurate surface temperature or climate sensitivity comes from three ingredients (Jeevanjee et al., 2022): (1) a realistic computation of CO₂, O₃, and water vapour radiative transfer; (2) the convective adjustment, where all unstable lapse rates are adjusted to the moist adiabatic lapse rate of −6.5 K m⁻¹; and (3) a fixed relative humidity in the radiative code to account for the water vapour feedback. Manabe and Strickler (1964) introduced the convective adjustment and Manabe and Wetherald (1967) the fixed relative humidity. Cloud microphysics was already present in the Manabe and Strickler (1964) model. Sarachik (1978) added a flux between the atmosphere and the ocean and found that in a model with no clouds and a surface with no heat capacity, the absence of the ocean compensates for the neglect of clouds. Rennó et al. (1994) explicitly computed the moisture profile with a cumulus convection scheme, leading to a wider model sensitivity to parameter changes. Recent developments in radiative–convective modelling focus on improving the cloud parameterization, creating a 3D radiative–convective model, or using it in a GCM (Wing et al., 2018). Radiative–convective models mainly follow the same tendency as GCM and ESM models. More and more complexity is added, and more and more parameters are present. All models are tuned with observations. Because of this tuning, we do not know whether these models are accurate in climates very far from the present Earth condition, such as palaeoclimates. The first “parameter” is the convection scheme in a radiative–convective model with a convective adjustment: the algorithm by which the lapse rate adjusts to the moist adiabatic can change the results. Also, small diffusion of temperature and moisture is used to stabilize the numerical scheme. Most importantly, the surface drag coefficient and wind speed are prescribed to values unknown in palaeoclimates, though they impact the surface energy budget. Indeed, the sensible and latent heat flux is usually calculated from gradients using a Fickian law, with a coefficient adjusted on observations. Since the physical (conservation) laws are the same for all climate models, the main differences between them are mostly linked to the tuning of these parameters: differences are therefore linked to different model settings, different choices of parametric formula, different tuning strategies, or different tuning data sets. It is probably one of the main problems in climate modelling.

To overcome this issue, we build a radiative–convective model with zero adjustable parameters in this study. Because the method used is completely new and has been little investigated, the model is built from scratch and appears to be a jump back into the 1970s. Many issues could be raised: the model is cloud-free (and we do not know how to compute clouds without parameters), it has a ground with no heat capacity and an infinite water reservoir, relative humidity is fixed to climatology in the radiative code and to 1 in the energy fluxes computation, albedo and solar constant are fixed, stationarity is assumed, and the model is 1D radiative–convective. However, the article's purpose is not to build a code as accurate as today's complex and highly tuned codes but to build a model as accurate as equivalent radiative–convective models built in the 1970s and 1980s. If it is possible to obtain similar or better results with a parameter-free model than a similar but tuned model, it would pave the way for creating a full climate model, as has been done for tuned models since the 1970s but this time not using tuning.

To get rid of parameters, the unknown variables are determined with a variational problem, the maximum of entropy production (MEP). The idea is to express the variational problem (entropy production and constraints) as a function of the unknown variables, like energy fluxes or water vapour fluxes, so they will adjust themselves to maximize entropy production. Therefore, we do not parameterize them. One could argue that MEP is just another method of parameterization. However, it is very different from the usual data assimilation techniques (see, for example, Lopez, 2007, for precipitation and clouds) because it does not use any observational data about the variable of interest to predict it (like temperature or precipitation). A MEP model was first used for climate by Paltridge (1975) to predict meridional fluxes and showed good agreement with observations. Unfortunately, it also had some parameterizations.

MEP is only a hypothesis and lacks rigorous mathematical proof. However, it seems very general and is used in many domains, like crystal growth, electric charge transfer, and biological evolution (Martyushev and Seleznev, 2006). First, the MEP hypothesis must not be mistaken for the second law of thermodynamics, which only states that entropy production is positive (σ≥0) but not necessarily maximized. Next, when equilibrium thermodynamics and entropy are reinterpreted with the formalism of information theory, it is possible to obtain the main results of equilibrium thermodynamics (Jaynes, 1957) using the maximum entropy principle (MaxEnt). A possible way to understand the maximum entropy production hypothesis is to see it as the non-equilibrium or time derivative “equivalent” of MaxEnt. Following this idea, Dewar (2003, 2005) tried to prove the MEP hypothesis using MaxEnt. However, to quote Martyushev (2021), “Dewar's argument not only involves a number of non-obvious fundamental assumptions but also is non-rigorous and erroneous in a number of points”. The MEP hypothesis has also been related to other variational approaches in fluid dynamics or climate (Ozawa et al., 2003), like the maximum dissipation rate when temperatures are fixed (Malkus, 1956) or the maximum generation of available potential energy when in a steady state (Lorenz, 1960). We do not try to demonstrate the MEP hypothesis here, but we prefer modelling the climate and the water cycle using it. Suppose results happen to be close to observations. In that case, the explanation as to why it works in our particular case will be left as an open question for theoreticians (see Martyushev, 2021, for a recent general review of the MEP hypothesis).

In climate, it has been used to predict oceanic or atmospheric horizontal heat fluxes (Grassl, 1981; Gerard et al., 1990; Lorenz et al., 2001; Paltridge et al., 2007; Herbert et al., 2011), as well as vertical heat fluxes (Ozawa and Ohmura, 1997; Pujol and Fort, 2002; Wang et al., 2008; Herbert et al., 2013) or both horizontal and vertical (Pascale et al., 2012), where always the entropy production of heat transfers due to atmospheric turbulence is maximized under some constraints. A limitation of using the MEP hypothesis (and variational formulations in general) is that all the variables of the variational problem need to be solved at once, making it difficult to add new phenomena. Another limitation is that to get meaningful results, we have to put on the variational problem constraints that are physically relevant and represent the main processes of the atmosphere (Goody, 2007; Dewar, 2009). This requires a lot of thinking and talent and is only sometimes easily solvable. Moreover, as the driver for heat transfers is the input of radiative energy in the system, the radiative code must be accurate to have results close to observations. Otherwise, we are limited to qualitative (Lorenz et al., 2001) but not quantitative (Goody, 2007) agreement with data. In previous studies, MEP has generally been used for straightforward cases or with additional parameters. For 1D vertical models, a grey atmosphere was used, and gravity was not considered. Recently, a new radiative–convective model was created with a more realistic radiative code (Herbert et al., 2013). Moreover, geopotential heat and latent heat were added (Labarre et al., 2019). This model already produces reasonable temperatures. This model is not based on a convective adjustment like Manabe and Wetherald (1967) but has the same crucial physical ingredients, allowing it to give similar results. However, the dynamical part (i.e. non-radiative) is treated without any adjustable parameter, giving the MEP model more universality because no one has ever told the model which values to take in any standard configuration. If the results are close to observations, it alerts us that there might be something powerful in the MEP hypothesis or the imposed constraints. In this study, the same model was used, and a new constraint on the water cycle was added, leading to a prediction of precipitation. The fact that this procedure can produce suitable water vapour fluxes has never been shown before.

2 The radiative code

The radiative code used here is the one of Herbert et al. (2013) (see their supplemental material for details) and is based on the net exchange formulation (Dufresne et al., 2005). It is more advanced than the grey atmosphere models used in previous studies (Ozawa and Ohmura, 1997), thus leading to results closer to observations. Let us consider an atmosphere divided into n+1 layers on the vertical axis (see Fig. 1). The net radiative energy budget ℛ_i (i.e. the input of energy thanks to radiation) is written as

$$\begin{array}{}\text{(1)}& {\mathcal{R}}_{\mathrm{i}}(T,q,{\mathrm{O}}_{\mathrm{3}},{\mathrm{CO}}_{\mathrm{2}},\mathit{\alpha })={\text{SW}}_{\mathrm{i}}(q,{\mathrm{O}}_{\mathrm{3}},\mathit{\alpha })+{\text{LW}}_{\mathrm{i}}(T,q,{\mathrm{CO}}_{\mathrm{2}}),\end{array}$$

where SW_i and LW_i are the solar and infrared net energy budgets. $q=M_{\mathrm{water}}/M_{\mathrm{air}}$, O₃, and CO₂ are prescribed vertical profiles of specific humidity and ozone and carbon dioxide concentrations, given on page 3 of McClatchey (1972), corresponding to a standard atmosphere, based on averaging observations. To take into account the water vapour feedback with temperature, the relative humidity profile $h=q/q_{\mathrm{s}}\left(T\right)$ is fixed for the computation of ℛ_i; therefore q varies with T. α is the albedo of the surface. T is the temperature profile. The parameters h, O₃, CO₂, and α are fixed, so the energy budget ℛ_i is a function of the temperatures only. Note that $T=\mathit{\{}{T}_j{\mathit{\}}}_{j=\mathrm{0},\mathrm{\dots },n}$, where T_j is the temperature in box j; hence ℛ_i(T) is a functional of the temperatures, i.e. a non-local function of T. That will be important for the variational problem.

Usually, when computing the radiative energy budget, one uses a local description of the radiative energy fluxes. In the net exchange formulation (NEF) framework (Dufresne et al., 2005), the description of energy transfer is global. Each radiative energy input in the layer i is broken down into the individual contributions of all atmospheric layers. The net energy exchange rate between layer i and j is ψ_ij. At this point, this formulation is strictly equivalent to the usual one. However, it makes it easier to develop approximations that reduce computational time while automatically satisfying the fundamental laws of physics: the energy exchange rate is antisymmetric (${\mathit{\psi }}_{ij}=-{\mathit{\psi }}_{ji}$). Thus, the total energy is conserved (${\sum }_{i,j}{\mathit{\psi }}_{ij}=\mathrm{0}$), and the second law of thermodynamics is also satisfied by keeping ψ_ij the same sign as T_j−T_i. Yet, because the resolution of the variational problem is strongly sensitive to the constraints imposed, it is of utmost importance that the laws of physics are rigorously satisfied.

To approximate the radiative transfer, the infrared spectrum is divided into 22 narrow bands, and the absorption coefficient for water vapour and carbon dioxide is calculated with the Goody (1952) statistical model and the data from Rodgers and Walshaw (1966). The diffusive approximation is made with the diffusion factor μ=1.66 for spatial integration. In the visible domain, water vapour and ozone absorption are computed with the parameterization of Lacis and Hansen (1974).

3 The maximization of entropy production

3.1 Previous setting: energy conservation

This constraint was used alone in Herbert et al. (2013). The atmosphere is divided into n+1 layers on the vertical axis. Each layer i has a temperature T_i, a variable of the variational problem. Between layer i and i+1, there is a non-radiative energy flux F_i, whose nature is not explicit yet (with no additional constraints, it may be interpreted as conduction). When considering a steady state, the total energy balance reads

$$\begin{array}{}\text{(2)}& {\mathcal{R}}_{\mathrm{i}}=F_{i+\mathrm{1}}-F_{\mathrm{i}},\end{array}$$

where F₀ and F_n+1 are taken equal to zero as if no energy (other than radiative) goes to the space or comes from the ground. The entropy production associated with these energy fluxes is written as

$$\begin{array}{}\text{(3)}& \mathit{\sigma }=-\sum _{i=\mathrm{1}}^n{F}_i\left({\displaystyle \frac{\mathrm{1}}{T_{i-\mathrm{1}}}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{i}}}}\right)\stackrel{\text{Eq.}\left(\mathrm{2}\right)}{=}\sum _{i=\mathrm{0}}^n-{\displaystyle \frac{{\mathcal{R}}_{\mathrm{i}}\left(T\right)}{T_{\mathrm{i}}}}.\end{array}$$

At stationary state, the total input of radiative energy must be equal to zero:

$$\begin{array}{}\text{(4)}& \sum _{i=\mathrm{0}}^n{\mathcal{R}}_{\mathrm{i}}\left(T\right)=\mathrm{0}.\end{array}$$

The entropy production (Eq. 3) is maximized under the constraint of energy conservation (Eq. 4), leading to the following variational problem:

$$\begin{array}{}\text{(5)}& \underset{(T_{\mathrm{0}},\mathrm{\dots },T_n)}{\text{max}}\left\{\sum _{i=\mathrm{0}}^n-{\displaystyle \frac{{\mathcal{R}}_{\mathrm{i}}\left(T\right)}{T_{\mathrm{i}}}}\left|\sum _{i=\mathrm{0}}^n{\mathcal{R}}_{\mathrm{i}}\left(T\right)=\mathrm{0}\right.\right\}.\end{array}$$

It is exactly the same problem as Eq. (24) of Herbert et al. (2013). The predicted lapse rate then has the correct magnitude in the lower troposphere but not in the upper layers (see blue square points in Fig. 2a).

Figure 1 Scheme of the 1D vertical box model, with the box 0 being infinitely small and other boxes separated by equal level of pressure. ℛ_i is the radiative energy balance and depends on all temperatures. F_i is the energy flux between adjacent boxes. m_i is the air mass flux. e_i is the specific energy (and depends on the local temperature), T_i is the temperature, and $q_{{\mathrm{S}}_{\mathrm{i}}}$ is the relative humidity taken at saturation (depends on the local temperature).

3.2 Previous setting: air convection

Air convection was studied in Labarre et al. (2019). So far in the article, the atmosphere's internal energy or air mass fluxes have not been expressed. It might give more physical results to put a constraint on how air masses transport the energy in the atmosphere in the variational problem. By taking into account the sensible heat, the gravitational energy, and the latent heat, the specific energy (the energy per unit mass) is written as

$$\begin{array}{}\text{(6)}& e_{\mathrm{i}}=C_p{T}_{\mathrm{i}}+g{z}_i+L{q}_i,\end{array}$$

where C_p is the heat capacity of the air, g the standard acceleration due to gravity, z_i the elevation of layer i, L the latent heat of vaporization, and q_i the specific humidity of water vapour ($m_{\mathrm{water}}/m_{\mathrm{air}}$). The elevation z_i is expressed as a function of the temperatures below (see Appendix A of Labarre et al., 2019), and the specific humidity q_i is taken equal to its value at saturation, q_i=q_s(T). Consequently, the specific energy e_i(T) is, like ℛ_i(T), a functional of temperature.

Now, take a mass flow rate m_i between layer i−1 and i. Suppose the air is transported adiabatically and then thermalizes once in layer i. It means that the energy in the air mass taken from layer i−1 is entirely transported to layer i. Thus, the upward energy transported equals m_ie_i−1. To easily conserve the air in all boxes, the same amount of air m_i is adiabatically taken from layer i to layer i−1, transporting downward energy equal to m_ie_i. Thereby, the total energy flux between layer i−1 and i is written as

$$\begin{array}{}\text{(7)}& F_{\mathrm{i}}=m_{\mathrm{i}}(e_{i-\mathrm{1}}-e_{\mathrm{i}}).\end{array}$$

Of course, for the reasoning to be consistent, we must have m_i≥0, which gives a new constraint on the energy fluxes and the temperatures. The following equation summarizes the variational problem:

$$\begin{array}{}\text{(8)}& \underset{(T_{\mathrm{0}},\mathrm{\dots },T_n),(m_{\mathrm{1}},\mathrm{\dots },m_n)}{\text{max}}\left\{\begin{array}{l}\sum _{i=\mathrm{0}}^n-\frac{{\mathcal{R}}_{\mathrm{i}}\left(T\right)}{T_{\mathrm{i}}}\left|\sum _{i=\mathrm{0}}^n{\mathcal{R}}_{\mathrm{i}}\right(T)=\mathrm{0}\\ \text{ and }{m}_{\mathrm{i}}\ge \mathrm{0}\\ \\ \text{with }{\mathcal{R}}_{\mathrm{i}}=F_{i+\mathrm{1}}-F_{\mathrm{i}},\\ F_{\mathrm{i}}=m_{\mathrm{i}}\left(e_{i-\mathrm{1}}\right(T)-e_{\mathrm{i}}(T\left)\right)\end{array}\right\}.\end{array}$$

It is the same variational problem as Eq. (11) of Labarre et al. (2019) (although expressed as a function of temperature instead of energy fluxes). The additional constraint allows for a much better simulation of the upper atmosphere (see orange diamonds Fig. 2).

3.3 New setting: water transport and precipitation

The model's novelty described here comes from the addition of water transport and precipitation. In the formulation above, water vapour is a function of the temperature only and thus has no reason to be conserved when transported by the air masses. Infinite levels of water vapour could be created or disappear. In this study, we add a constraint on conserving water vapour that is supposed to mimic precipitation. We impose that water vapour cannot appear when transported, but it can disappear, and we call this phenomenon precipitation as if water vapour was transformed into liquid water. The flow rate m_i transports upward between layer i−1 and i an amount of water equal to m_iq_i−1, where q_i is the specific humidity of water vapour in the air, and it transports downward an amount of water equal to m_iq_i. Then the water flux between layers i−1 and i is written as (similarly to Eq. 7)

$$\begin{array}{}\text{(9)}& F_{\mathrm{i}}^{\mathrm{w}}=m_{\mathrm{i}}(q_{i-\mathrm{1}}-q_i).\end{array}$$

The amount of water vapour that disappears in layer i is written as

$$\begin{array}{}\text{(10)}& P_{\mathrm{i}}=F_{i+\mathrm{1}}^{\mathrm{w}}-F_{\mathrm{i}}^{\mathrm{w}},\end{array}$$

where $F_{\mathrm{0}}^{\mathrm{w}}$ and $F_{n+\mathrm{1}}^{\mathrm{w}}$ are taken equal to zero as if no water comes from underground or goes to space. The layer i=0 is a surface boundary layer, supposed to be very thin, and plays the role of the surface. For $i=\mathrm{1},\mathrm{\dots },n$, we impose that P_i≥0, and P_i is called precipitation. On the layer i=0, because ${\sum }_{\mathrm{i}}{P}_{\mathrm{i}}=\mathrm{0}$ we have

$$\begin{array}{}\text{(11)}& P_{\mathrm{0}}=-\sum _{i=\mathrm{1}}^n{P}_{\mathrm{i}},\end{array}$$

where −P₀ is the evaporation in layer i=0 and is equal to the total precipitation. The specific humidity q_i is considered equal to its value at saturation q_s(T) and then depends only on the temperature. The variational problem can be summarized by

$$\begin{array}{}\text{(12)}& \underset{(T_{\mathrm{0}},\mathrm{\dots },T_n),(m_{\mathrm{1}},\mathrm{\dots },m_n)}{\text{max}}\left\{\begin{array}{l}\sum _{i=\mathrm{0}}^n-\frac{{\mathcal{R}}_{\mathrm{i}}\left(T\right)}{T_{\mathrm{i}}}\left|\sum _{i=\mathrm{0}}^n{\mathcal{R}}_{\mathrm{i}}\right(T)=\mathrm{0}\\ \text{ and }{m}_{\mathrm{i}}\ge \mathrm{0}\text{, }{P}_{\mathrm{i}}\ge \mathrm{0}\text{ for }i\ge \mathrm{1}\\ \\ \text{with }{\mathcal{R}}_{\mathrm{i}}=F_{i+\mathrm{1}}-F_{\mathrm{i}}\text{, }{F}_{\mathrm{i}}\\ =m_{\mathrm{i}}\left(e_{i-\mathrm{1}}\right(T)-e_{\mathrm{i}}(T\left)\right)\text{, }\\ P_{\mathrm{i}}=F_{i+\mathrm{1}}^{\mathrm{w}}-F_{\mathrm{i}}^{\mathrm{w}}\text{, }{F}_{\mathrm{i}}^{\mathrm{w}}\\ =m_{\mathrm{i}}\left(q_{S_{i-\mathrm{1}}}\right(T)-q_{S_{\mathrm{i}}}(T\left)\right)\end{array}\right\}.\end{array}$$

Although it looks like a complex equation, it is, in fact, very simple when one considers that absolutely all the physics of the model is contained in it.

4 Numerical resolution

The variational problems (Eqs. 5, 8, and 12) are solved using a sequential quadratic programming algorithm (Kraft, 1988; Virtanen et al., 2020). The basic principle is that it takes some initial conditions and performs a gradient descent until it finds a local maximum. Such an algorithm mathematically converges to a global maximum for a convex problem; however, for a non-convex problem, there is no guarantee that the maximum found is global. Several techniques that are not detailed here are used to get better results. One possibility is to manually test different initial conditions and see which gives the highest entropy production. For a given problem, every result presented here is the one with the highest entropy production found. Given that we tested many different initial conditions and found only a few local maxima, we are confident our results represent a global maximum. Nevertheless, it is still possible that another better result mathematically exists.

5 Results

We solve the problems with 21 boxes (1 surface boundary layer with albedo α and 20 atmospheric boxes; see Fig. 1), with prescribed standard vertical profiles of O₃ and CO₂ concentrations and relative humidity taken from McClatchey (1972), corresponding to a tropical atmosphere. Indeed, because it is a 1D vertical model, taking midlatitude or boreal profiles would not make much sense, as horizontal fluxes become significant there. In the radiative code, the relative humidity $h=q/q_{\mathrm{S}}\left(T\right)$ is fixed; therefore more H₂O is present when temperature increases, which is an important positive feedback for climate sensitivity. However, in the variational problem, q_i is still equal to q_s(T) (Eq. 6). This decoupling may be inconsistent, but it is common in the climate community or GCMs to examine models decoupled before studying fully coupled model configurations. We also choose to proceed in this way to keep the problem simpler and to be able to interpret results more easily. As in Rennó et al. (1994), incoming solar radiation is fixed to 342 W m⁻², and albedo is equal to α=0.1.

Results for Eqs. (5), (8), and (12) are shown Fig. 2. For temperature and energy, the lowest point corresponds to the surface. For comparison, standard temperatures for a tropical atmosphere (based on observations) are plotted in Fig. 2a, taken from McClatchey (1972), as well as mean temperatures from the IPSL-CM6A-LR model between 23° S and 23° N. The observations are not taken in a specific area of the tropics. However, if we knew a place matching the requirements of the model, that is to say, where horizontal fluxes are negligible, and where we know the greenhouse gas concentration, it could be used and compared to observations at this place. The comparison is made here for spatially average tropical conditions, using standard greenhouse gas concentrations and temperatures, as is customary for radiative–convective models. Also, we only consider a stationary atmosphere, so the results should be interpreted only in a climatological context, not a meteorological one. The comparison should remain qualitative, as our model's purpose is not to fit the observations precisely but to give a “proof of concept” that obtaining relevant results with the MEP hypothesis is possible. First, when putting more and more constraints into the variational problem, it is expected that the entropy production found decreases because the space of possibilities wanes. We obtain the following: the entropy production is equal to σ=53.9 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ for Eq. (5) where only energy is conserved, σ=44.3 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ for Eq. (8) where a specific pattern of mass fluxes transports energy, and σ=41.1 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ for Eq. (12) where water vapour is not allowed to appear.

5.1 With only energy conservation (blue square points)

The use of energy conservation only has already been studied in Herbert et al. (2013). Although the numerical resolution has been improved, the results here are very similar to Herbert et al. (2013), but we show them for comparison with the next ones. The only variables present in Eq. (5) are temperatures T and energy fluxes F. However, afterwards it is possible to compute variables like energy e (Eq. 6), mass fluxes m (Eq. 7), or precipitation P (Eq. 10), but of course, there is no reason for m or P to obey the constraint of positivity, and they do not. When looking at the entropy production equation (Eq. 3), one sees two terms: an energy flux and a gradient of inverse temperature. To maximize the entropy production, variables find a balance between a state where energy fluxes are very high but temperatures homogeneous (so σ=0) and a state where temperature gradients are high but energy fluxes are equal to zero (also σ=0). Between the two, a balance is found where entropy production is maximized. Results are plotted in Fig. 2a for temperature and Fig. 2b for energy fluxes. As expected, temperatures are higher near the ground because the atmosphere is essentially heated by the solar radiation intercepted at z=0. The energy fluxes compensate a bit for the temperature difference by going upwards (i.e. are positive Fig. 2b). The specific energy (Eq. 6, Fig. 2c) decreases at the bottom because sensible heat C_pT is the dominant term and increases at the top because gravity gz becomes the dominant term. “Mass fluxes” are computed with Eq. (7). Because F is always positive, the sign of m is always opposite to that of the energy gradients ∇e. When these gradients are negative, m is positive, but higher in the atmosphere when energy gradients become positive, the computed m is then negative, which has no physical meaning regarding Eq. (7). Then, imposing m>0 becomes natural, as done in the following paragraphs.

5.2 With a pattern of convection (orange diamond points)

The use of energy and mass conservation has already been studied in Labarre et al. (2019). Because the numerical resolution has been improved, results here differ slightly from Labarre et al. (2019). Adding the constraint of convection (Eq. 7) forces the energy fluxes (Fig. 2b) to be opposed to the energy gradients. Thus, energy (Fig. 2c) gradients must remain negative to keep positive upward energy fluxes (energy must decrease with height). To do so, temperatures (Fig. 2a) adapt themselves to counteract the gravity, leading to a zone in the middle of the atmosphere where energy gradients equal zero and mass fluxes (Fig. 2d) equal +∞; the atmosphere is perfectly mixed in this area, energy follows the moist adiabatic lapse rate (Fig. 2c, a vertical line). Thus, precipitation (Fig. 2e) is infinitely positive or negative. Higher in the sky, geopotential becomes too strong, energy gradients become positive, and energy fluxes equal 0. This region defines the stratosphere at about P≈250 hPa and z≈9 km, where no convection occurs. By adding a convection pattern and maximizing entropy production, we see a moist adiabatic lapse rate in the middle of the atmosphere, and a stratosphere naturally emerges. This model already gives suitable results for temperatures as the profile is similar to the IPSL-CM6A-LR model (light blue triangles) or the McClatchey (1972) standard tropical atmosphere (circles). The temperature profile remains close when adding a new constraint to compute precipitation.

Figure 2 Three different cases: (1) Eq. (5), (2) Eq. (8), and (3) Eq. (12). (1) σ=53.917, (2) 44.304, and (3) 41.108 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$.

5.3 With a constraint on precipitation (grey triangle points)

Adding the constraint on precipitation (Eq. 10) is the article's novelty. Nevertheless, the philosophy remains the same: because the atmosphere is heated from below, energy fluxes (Fig. 2b) want to go upward, but they have to be opposed to the energy gradients (Fig. 2c), so these gradients want to be negative or to go to zero, leading to infinite mass fluxes (Fig. 2d) in the middle of the atmosphere. However, saying that no water vapour can be created prevents infinite mass fluxes and nil energy gradients. Consequently, energy decreases (while energy fluxes are positive) below a zone we call the tropopause; above, energy increases (while energy fluxes are equal to zero). Regarding precipitation (Fig. 2e), some water vapour is taken in box number 0 and continues going up to maximize mass fluxes. Then it reaches the tropopause, where mass fluxes become nil because there is not enough energy to go up, and then water vapour disappears; i.e. it precipitates.

The computed precipitation is equal to 1.2 m yr⁻¹, which is the correct order of magnitude of tropical precipitation. Comparison of this result with real-world or modelled tropical precipitation depends on the box size chosen for the tropical area. For example, the average Earth system model IPSL-CM6-LR shows that precipitation between −23 and +23° of latitude is 1.4 m yr⁻¹. For real-world data, average precipitation between 30° S and 30° N between 1980 and 1994 is 1.3 m yr⁻¹ (Fig. 8 of Xie and Arkin, 1997), and zonally averaged precipitation between 1979 and 2001 lies between 0.6 and 2 m yr⁻¹ (Fig. 5 of Adler et al., 2003). Note that the MEP model gives less precipitation, and the other model gives more precipitation than observations. We will see that it is a common bias for models to overestimate precipitation.

Another local maximum of entropy production can give, for example, 2.2 m yr⁻¹, though with less overall entropy production. This “uncertainty” in computed precipitation is probably due to the low resolution of only 20 boxes because choosing one box or another to precipitate leads to a different value of precipitation. To test this hypothesis, it is possible to simplify the problem by imposing the precipitation value in only one box. Then, the problem is numerically easily solved, and precipitation converges when the number of boxes increases. With 81 boxes, the computed precipitation is 1.15 m yr⁻¹, and the other local maximum of entropy production leads to precipitation of 1.21 m yr⁻¹. Thus, the 1.2 m yr⁻¹ precipitation is a robust value.

To the author's knowledge, it is the first time precipitation is computed with a model using maximizing entropy production and without any data-tuned parameter. At least, the fact that the good order of magnitude of precipitation can be computed with little knowledge is of prime theoretical importance for climate scientists because it means the radiative transfer, or greenhouse gases, mainly drives atmospheric precipitation. This statement is not new and can also be deducted from Pierrehumbert (2002), Rennó et al. (1994), or Eq. (9) of O’Gorman and Schneider (2008). Still, variation in surface sensible heat flux or how the surface radiative transfer depends on other processes is a source of discordance between models, leading to slightly different computed surface latent heat flux (or precipitation). Our model has less bias than radiative–convective models or the one from the atmospheric model intercomparison project, and further investigation is needed to understand if it is by chance or if there is some profoundly hidden physics in the MEP hypothesis or the imposed constraints.

5.4 Doubling the CO₂ concentration

When looking at absolute values, the MEP approach seems to provide a good order of magnitude for climate variables. It is interesting to test if it can capture small changes in external forcing. A classic test for climate models is to look at the climate sensitivity, which is the difference of temperature at 1.50 m between conditions where CO₂ is at pre-industrial level (280 ppm) and conditions where CO₂ is doubled (560 ppm). The only feedbacks present here are the water vapour feedback incorporated by fixing the relative humidity in the radiative code, and the lapse rate feedback. Because box number 0 is a thin surface boundary layer, we take the temperature of layer number 1, whose middle is at 988 hPa (∼220 m). The climate sensitivity is 1.1 K for Eq. (5), 0.7 K for Eq. (8), and 1.0 K for Eq. (12) (see Fig. 3). For comparison, the climate sensitivity of the state-of-the-art IPSL-CM6A-LR model is 3.0 K, which is much more. The limited sensitivity to CO₂ is not surprising since our model only represents a few amplifying phenomena. For instance, the albedo is fixed, there is no deep atmospheric convection, and the ground has no heat capacity. The ground is an infinite water reservoir, the model is 1D, and no clouds are in the atmosphere. How to integrate (some of) these additional processes in the MEP procedure remains to be investigated, as well as how this would affect climate sensitivity.

Nevertheless, the model of Manabe and Wetherald (1967) contains physics similar to the MEP model and predicts 2.9 K (see Table 5). The result looks better because the convective adjustment method used in Manabe and Wetherald (1967) is very efficient in transferring a rise of temperature from the top of the atmosphere into the bottom of the atmosphere, and that is why the climate sensitivity, which only regards the temperature at the bottom of the atmosphere, is so well predicted. In the MEP case, the temperature at the bottom of the atmosphere is much less constrained by radiative transfer, which is why the result seems less good. However, looking not at one temperature but at temperature vertical profiles, the Manabe and Wetherald (1967) model computes a constant elevation of temperature in the whole troposphere (see their Fig. 16), which is not at all consistent with the vertical profile computed by the state-of-the-art IPSL-CM6A-LR model (SSP245-CTRL Fig. 3). Interestingly, radiative–convective models that do not use a mean moist convective adjustment (of 6.5 K m⁻¹) but a local moist adiabatic convective adjustment or a cumulus model (Lindzen et al., 1982) predict both a better temperature profile and a considerably reduced climate sensitivity. For example, depending on the imposed parameters, Lindzen et al. (1982) find a climate sensitivity of 1.4 and 1.7 K with a moist adiabatic convective adjustment and 1.6, 2.2, 2.3, or 0.9 K with a cumulus model (see Tables 1 and 2 of Lindzen et al., 1982).

Figure 3 Differences of temperatures profiles between CO₂ at 180 ppm (or 560 ppm) and CO₂ at 280 ppm.

The vertical temperature distribution within the atmosphere is plotted for Eqs. (5), (8), and (12) (Fig. 3), that is the difference of temperatures between a case with 560 ppm (or 180 ppm corresponding to the last glacial maximum) and a case with 280 ppm. The difference between 560 and 280 ppm is also plotted for the IPSL-CM6A-LR model (SSP245-CTRL), with a temperature average taken in the tropical region (between ±23°). When looking at the shape of temperature distribution and comparing it to the SSP245-CTRL, the best model seems to be the problem (Eq. 8) since the temperature difference increases with height and then decreases in the stratosphere. For the Eq. (12) model, it just decreases in the entire atmosphere. The less-constrained CONV experiment provides a better sensitivity of the temperature profile to CO₂ doubling and, in that sense, does better than the Eq. (12) case or Manabe and Wetherald (1967) model. Results from more evoluted radiative–convective models are similar to Eq. (8) (see, for example, Figs. 3 and 8 of Lindzen et al., 1982). Considering the simplicity of the model and the absence of parameterizations, it was expected that correct values for all model outputs would not be obtained. The fact that Eq. (8) is better and compares well with similar parameterized models suggests that the Eq. (12) experiment is over-constraining the water cycle, something that could be relaxed by introducing new degrees of freedom. For example, the convection pattern could allow mass fluxes between layers that are not adjacent, allowing deep convection to occur, or water vapour could be allowed to vary between zero and saturation, or the setup could be extended from 1D to 2D. We think that once this is done, the sensitivity of the temperature profile to CO₂ should get close to the Eq. (8) experiment. Unfortunately, this complexifies the optimization problem, and we have yet to be able to test this idea.

Moreover, when interpreting these results, one should remember that depending on the resolution method of the variational problem (or the choice of initial conditions), results may differ by about 1 K. They differ even more by choosing an arbitrary local maximum of entropy production instead of the “supposedly” global maximum. Indeed, we note that because we changed the resolution method and found a new result with higher entropy production, the climate sensitivity of the problem (Eq. 8) is a bit different than in Fig. 6 of Labarre et al. (2019). However, the purpose of this model was not to produce precise or reliable predictions. More modestly, this model aims to demonstrate that some relevant elements of the climate system (temperature and precipitation) can be computed with only a minimal set of hypotheses, i.e. conservation laws only, without any parameter tuning.

5.5 Fixing the relative humidity

In the above simulations, in the energy equation (Eq. 6), the specific humidity q has been fixed to its value at saturation q_s(T). In other words, the relative humidity is fixed to $h=q/q_{\mathrm{s}}\left(T\right)=\mathrm{1}$. Relaxing this constraint by letting the relative humidity vary between 0 and 1 poses a numerical problem that has yet to be solved. However, one may wonder how fixing the relative humidity by hand would change the results. In this section, we solve the problem with a constraint on precipitation (Eq. 12) and a relative humidity constant in the entire column but chosen between 0.1 and 1. h=1 corresponds to the Eq. (12) setting. As before, a possible method is to test different initial conditions and see which results give the maximum entropy production. However, this method seems to fail when h<0.5, and the results are not robust anymore. Fortunately, phenomenology shows that in all our simulations, the case where precipitation happens in only one box seems to always correspond to the maximum entropy production. With this knowledge, we use a simpler way to obtain the MEP solution. We use a method where the place and value of precipitation are imposed to obtain the result more directly. We find the maximum entropy production by exploring the phase space by changing the value of P. We use the following protocol: relative humidity (h), precipitation (P), and the box where precipitation occurs (n_p) are fixed. A value of entropy production σ is found for each set of $(h,P,n_{\mathrm{p}})$. The explored phase space is $h\in \mathit{\{}\mathrm{0.1},\mathrm{0.2},\mathrm{\dots },\mathrm{1}\mathit{\}}$, $P\in \mathit{\{}\mathrm{0.1},\mathrm{0.2},\mathrm{\dots },\mathrm{10}\mathit{\}}$, $n_{\mathrm{p}}\in \mathit{\{}\mathrm{13},\mathrm{14},\mathrm{\dots },\mathrm{19}\mathit{\}}$. The number of atmospheric boxes is fixed at n=20. The results are shown in Fig. 4, with the obtained values of maximum entropy production σ and precipitation P. For all values of h, the maximum of entropy production corresponds to precipitation in box number n_p=15, as in the Eq. (12) problem (corresponding to h=1 here).

Figure 4 Values of maximum entropy production σ and corresponding precipitation P for each constant profile of relative humidity h. The maximum only corresponds to precipitation in box number 15 (with n=20). The maximum between all h is h=1, with σ=41.108 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ and P=1.2 m yr⁻¹.

When h tends to zero, precipitation seems to go to infinity. We stopped at P=10 m yr⁻¹, so the true value of P is greater than 10 m yr⁻¹ for h≤0.3. Further analysis shows that when h tends to zero, the function σ(P) has less and less of a clear maximum and becomes very flat, so very different values of precipitation can give slightly different values of entropy production. This explains why results were not robust when using the classic convergence algorithm when h is low: the phase space becomes very flat in the region of the maximum. The important result arising from Fig. 4 is that the maximum of entropy production corresponds to h=1, and P=1.2 m yr⁻¹. This means that in this 1D vertical and stationary case, the system is the most efficient to create entropy when water vapour is saturated. We did not test all possible relative humidity profiles, but we expect that the MEP solution is close to h=1 in each box.

6 Discussion

The state-of-the-art GCMs and ESMs and MEP models are based on the same conservation laws. In a GCM, the local conservation laws lead to partial derivative equations that are true in the limit of infinitely small differentials. The momentum conservation is the Navier–Stokes equation (present only to the horizontal), the energy conservation is the thermal energy equation, and the mass conservation is $\mathrm{\nabla }\cdot u=\mathrm{0}$. These equations were demonstrated for infinitely small increments. However, in GCMs and ESMs, they are integrated under a grid that needs to be smaller. Indeed, because of non-linearity, small scales do have an impact on large scales. Therefore, something more than the first conservation equations is needed for consistent results. Additional equations involving tunable parameters are added and are sometimes called “closure equations”. In the MEP framework, the energy conservation is Eq. (4), and the mass conservation is immediately imposed by the convection pattern (Eq. 6) and constraint m≥0. The water conservation is imposed by P≥0 in Eq. (12). So, with MEP, the conservation laws are defined as constraints of an optimization problem, and unknown variables are resolved simultaneously to reach the maximum entropy production.

Everything else in our MEP model is similar to what is done in a GCM. The radiative code is based on integrating Planck's law on different wavelength bands corresponding to different constant extinction coefficients (see supplementary materials of Herbert et al., 2013). The air is considered an ideal gas, and the hypothesis of hydrostatic equilibrium is made ($\mathrm{d}p=-\mathit{\rho }g\mathrm{d}z$). Of course, well-known parameters like heat capacity C_p or the relation between q_s and T (see Appendix B of Labarre et al., 2019) are not variables of the optimization problem but just taken equal to well-established values used in GCMs and ESMs.

Still, there are many reasons why our MEP model could give different results than an ESM like the IPSL-CM6A-LR. Our MEP model has no continental surface and no energy flux between the ground and the atmosphere. Because there is no constraint on evaporation at the surface level, the ground can be seen as an infinite water reservoir, like an ocean. The ground heat flux is crucial on daily or seasonal timescales. However, it must equal zero in a stationary state without sub-surface fluxes (neglecting oceanic horizontal fluxes for ocean surfaces and geothermal fluxes on continental ones). The ground heat flux will be important if we add time to the problem to see a daily or seasonal cycle. However, this is a work beyond the scope of this article.

Also, no clouds (i.e. liquid water) are present in the air, although they are known to have an important impact on the radiative forcing. However, according to Sarachik (1978), neglecting the clouds introduces a false cooling of the surface that compensates for the false cooling by neglecting the heat transfer in the ocean. Moreover, the model is only 1D vertical and works well only for a tropical column where vertical fluxes dominate. However, O’Gorman and Schneider (2008) showed that the radiative–convective approximation is good for computing average precipitation. To refine the results, it would be easy to use fixed prescribed horizontal fluxes, but this differs from our goal. Also, there is no theoretical difficulty in extending this model to 2D or 3D and obtaining horizontal fluxes, but it is currently very technically challenging.

Finally, a reason for getting different results could be the possible lack of validity of the MEP hypothesis. That said, obtaining precipitation values close to the IPSL-CM6A-LR model or observations by Xie and Arkin (1997), 1.2 m yr⁻¹ compared to 1.4 and 1.3 m yr⁻¹, is very encouraging. Intuitively, given the fact that the model is built in a fully non-parametric fashion and is so minimalistic, we could have expected much more discrepancy. However, O’Gorman and Schneider (2008) showed that radiative transfer strongly constrains precipitation (see their Fig. 3). For example, the world's mean precipitation cannot exceed 2.1 m yr⁻¹ (Eq. 6 of O’Gorman and Schneider, 2008). Indeed, there is a balance at the surface between the net radiative energy budget, sensible heat flux, and latent heat flux. In the tropics, the sensible heat flux is 1 order of magnitude less than the latent heat flux and can be neglected. The simplified energy budget allows us to give a precipitation value, given a temperature; see Fig. 3 of O’Gorman and Schneider (2008). If the surface temperature is well constrained, which is the case in a tuned parameterized model, precipitation is also well constrained. The MEP model reproduces this simplicity: both a dominant latent heat flux and a good surface temperature. At 300 K, O’Gorman and Schneider (2008) find precipitation of about 1.8 m yr⁻¹. The difference with us comes from the grey atmosphere radiative code they use that is too simple. Interestingly, they prove that the same results regarding precipitation obtained from a GCM can be obtained with a radiative–convective model (with the surface wind speed parameter chosen from the mean value of the GCM). This means the radiative–convective approach is a reasonable first approximation to investigate precipitation.

One could wonder if there is a need to add a constraint on water vapour conservation to compute precipitation. Indeed, Eqs. (5), (8), and (12) have similar surface net radiative energy budget (84, 82, and 98 W m⁻²), and saying sensible heat flux is negligible, one could already deduce the surface latent heat flux thus evaporation (i.e. precipitation). However, in the three models, the computed sensible heat flux is not always negligible: for Eq. (5), the surface sensible heat flux is 19 W m⁻², and the surface latent heat flux is 65 W m⁻². For Eq. (8), it is 15 and 67 W m⁻². And for Eq. (12), it is 2 and 96 W m⁻². Evaporation is higher in the Eq. (12) model. Most importantly, Eqs. (5) and (8) do not conserve water (precipitation defined as the convergence of water fluxes is either negative or infinite): evaluating the water fluxes not at the surface but at another model level would give completely arbitrary results.

The MEP model could be improved by exploring several approaches. First, the specific humidity of water vapour q could be chosen not equal to saturation. Then, it needs to be clarified if an additional constraint on precipitation should be imposed, for example, saying that precipitation can occur only if q=q_s, that is $(q-q_{\mathrm{s}})P=\mathrm{0}$, where P is precipitation. Moreover, this constraint is highly non-convex and numerically very difficult to solve. Second, convection is not allowed between every layer because air masses are compelled to move to adjacent layers. However, in the tropics, there is a phenomenon called deep convection, where an air mass can go adiabatically from the bottom to the top of the troposphere. This phenomenon is not authorized with the convection pattern imposed here. However, this could be implemented by changing Eq. (6) and adding fluxes for non-adjacent layers.

7 Conclusions

Since Ozawa and Ohmura (1997), many improvements have been made. Taking a more realistic radiative code (Herbert et al., 2013) leads to a stable atmosphere: the potential temperature increases with altitude. Adding a convection pattern (Labarre et al., 2019) gives much more realistic results in terms of temperature and reveals a stratosphere up to ≈250 hPa where no convection occurs. Then, imposing a constraint on water vapour conservation leads to precipitation as close as or closer to observations than a GCM or an ESM would find. These results look great in absolute value, but results seem less satisfying when performing a sensitivity experiment such as looking at the temperature difference when doubling CO₂. We hope that future work that adds degrees of freedom in the model should solve this issue. Nonetheless, for the first time, this article demonstrates that it is possible to compute tropical precipitation in a radiative–convective model only from conservation laws, replacing the usual atmospheric parameterizations with an optimization procedure (MEP hypothesis).

Several approaches, such as changing the convection pattern and letting relative humidity vary, need to be explored in the future. Along these lines, it might be possible to build a 2D or 3D climate model, with a representation of clouds, vegetation, and oceans, whose “closure equations” would not at all be based on parameters that are “tuned towards observations”. The way clouds could be added without using any parameter is complicated. However, letting the relative humidity vary between zero and saturation, adding fluxes between non-adjacent boxes (i.e. deep convection), or making the problem 2D or 3D is mathematically straightforward. However, in such more complex settings, we have to solve a non-linear and non-convex optimization problem with many more variables. Currently, technical difficulties appear in the algorithm's convergence: the less convex the problem is, the more variables and non-linearities there are, and the harder it is to find the absolute global maximum. This technical problem might be solved using more specific optimization algorithms or finding a more straightforward mathematical but equivalent formulation.

Such a model could be used in contexts where tuning is impossible, such as other planets or palaeoclimates. The only mandatory knowledge is the atmosphere's chemical composition, the value of the solar radiation, and the ground albedo. This simplicity opens new and exciting perspectives in palaeoclimate or exoplanet climate modelling.

Appendix A: A global and a local maximum

Figure A1 Two different cases: (1) a global maximum of entropy production σ=41.108 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ and (2) a local maximum of entropy production σ=40.078 $\mathrm{mW}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$.