A radiative-convective model based on constrained maximum entropy production

The representation of atmospheric convection induced by radiative forcing is a long-standing question mainly because turbulence plays a key role in the transport of energy as sensible heat, geopotential, and latent heat. Recent works have tried using the maximum entropy production (MEP) conjecture as a closure hypothesis in 1-D simple climate models to compute implicitly temperatures and the vertical energy flux. However, these models fail to reproduce realistic profiles. To solve the problem, we describe the energy fluxes as a product of a positive mass mixing coefficient with the corresponding energy gradient. This appears as a constraint which imposes the direction and/or limits the amplitude of the energy fluxes. It leads to a different MEP steady state which naturally depends on the considered energy terms in the model. Accounting for this additional constraint improves the results. Temperature and energy flux are closer to observations, and we reproduce stratification when we consider the geopotential. Variations in the atmospheric composition, such as a doubling of the carbon dioxide concentration, are also investigated.

1 Introduction

The climate system is complex and usually divided into different components: atmosphere, ocean, cryosphere, lithosphere, and biosphere (Peixoto and Oort, 1992). There are different approaches to climate modelling (Randall et al., 2007). We can classify them in three main classes. Global climate models (GCMs) are the more sophisticated ones (see Dufresne et al., 2013, for an example). They explicitly represent the circulation of the atmosphere and ocean. Earth models of intermediate complexity (EMICs) simulate the Earth system with more simplifications than GCMs (see Goosse et al., 2010, for an example). These simplifications allow simulations over larger time periods, which is useful to study past climates. Simple climate models (SCMs) use only a few key processes to answer specific questions (see Paillard, 1998, for an example).

Both the complex and simple models have different strengths and weaknesses and are used for different applications. For example, GCMs are largely used to make climate projections for the next century. Since the numerical resolution of dynamical equations from the micro-scale (of the order $\simeq {\mathrm{10}}^{-\mathrm{3}}$ m for viscous dissipation) to the scale of interest (≃10⁷ m for the typical size of the Earth) is still impossible, GCMs need to represent sub-grid processes such as small eddies, convection, or cloud formation. To do it, models usually express the intensity of fluxes due to unresolved phenomena as a function of the resolved variables. This approach is called a “turbulent closure” and usually requires the introduction of empirical parameters such as turbulent master length scale, turbulent velocity diffusion terms (Mellor and Yamada, 1982), or the use of different quantities like convective available potential energy (CAPE) or convective inhibition (CIN) to fix the convective intensity (Yano et al., 2013). These parameterizations change from one model to another, resulting in different predictions (Stevens and Bony, 2013). They also require adjusting the numerous free parameters (“tuning”) in order to track observations (Hourdin et al., 2017). SCMs appear as an interesting alternative when considering past or future climate over a larger period like glacial–interglacial cycles. Indeed, an SCM is a set of a reasonable number of equations and physical quantities, providing an easier assessment of the impact of the parameters of the models (such as the concentration of greenhouse gases for atmospheric models). Many SCMs are based on the idea that the computation of all the microscopic details may be unnecessary if we are interested in quantities at larger spatio-temporal scales.

Consequently, a lot of SCMs describe the Earth simply with energetic considerations (North et al., 1981). Those models are called energy balance models (EBMs). The atmosphere is mainly driven by radiative forcing: solar radiation imparts energy to the Earth, which emits infra-red radiation to space. This heating is not homogeneous around the globe for different reasons such as geometry or variations in albedo with the nature of the ground. The insolation is more important for the tropics than the poles and leads to latitudinal heat transport. For the vertical axis, the amount of radiative energy absorbed by the Earth naturally depends on the atmospheric components. The ground usually receives more solar radiation because of the relative transparency of the atmosphere. As a result, the atmosphere is heated from below, which may lead to unstable situations where the temperature gradient exceeds the adiabatic gradient. Then, this causes atmospheric motion and vertical heat transport named convection. Since EBMs are usually based only on the energy budget, it is necessary to specify a relation between the energy fluxes and the temperature gradient, called a closure hypothesis. So EBMs mainly differ by the representation of the fluxes between the fluid layers of the Earth (atmosphere and ocean). For example, horizontal fluxes are sometimes represented by purely diffusion terms (North et al., 1981). On the other hand, vertical energy transport has been modelled using different approaches such as the convective adjustment. The latter consists of computing the temperature profile at radiative equilibrium for stable regions and adjusting it where the critical gradient is exceeded (Manabe and Strickler, 1964). Representing both the horizontal and the vertical energy fluxes is an important issue for EBMs. This concerns the subject of the present paper.

Since the 70s (Paltridge, 1975), maximum entropy production (MEP; Martyushev and Seleznev, 2006) is also used as a closure hypothesis in EBMs. This conjecture stipulates that the climatic system (or one of its component) optimizes its entropy production due to internal heat transfers. It allows for computing implicitly (i.e. without the computation of the dynamics) horizontal fluxes (O'brien and Stephens, 1995; Lorenz et al., 2001) and vertical fluxes (Ozawa and Ohmura, 1997; Pujol and Fort, 2002) without the parameterizations required in more conventional models. Former MEP-based models (MEPMs) have been criticized for three main reasons. One is the absence of dynamics and the validity of MEP (Rodgers, 1976). The second criticism deals with the extra parameterizations or the assumptions used in MEPMs. Indeed, one may ask oneself if the successes of the models are really due to the MEP hypothesis or to tuning or other ingredients (Goody, 2007). The final criticism concerns the usually simplified description of the radiative forcing in these models. Recently, a MEPM overcoming the last two criticisms has been built in Herbert (2012). It includes a refined description of the radiative budget in the net exchange formalism, without extra assumptions. The only adjustable quantities concern the radiative budget, such as the albedo, and not the atmospheric or oceanic energy transport. The model provides a relatively good approximation for the temperature and horizontal heat fluxes (Herbert et al., 2011b).

However, the vertical energy fluxes are still overestimated in such models when compared to observations or conventional radiative-convective models (RCMs) like Manabe and Strickler (1964). Furthermore, the energy fluxes are not always oriented against the energy gradient and it does not predict stratification in the upper atmosphere. This is not surprising because geopotential was not taken into account. Yet, we know from fluid mechanics that gravity plays a major role in natural convection (Rieutord, 2015). Gravity is also obviously responsible for stratification in the upper atmosphere. In this paper, we develop a MEPM that describes more properly the atmospheric convection. In the same spirit as the previous SCMs, we do not attempt to resolve the dynamical equations, but we add only some key features. Two ingredients are introduced to represent vertical heat fluxes more correctly. The first one is to describe energy transport as the product of a non-homogeneous mixing mass coefficient and the specific energy gradient. This brings a new constraint into a MEPM. The second one is to consider different energy terms: sensible heat, geopotential, and latent heat. We show that this simplified description of the energy transport, combined with the MEP closure hypothesis, can lead to relatively realistic results.

The outline of this paper is as follows. In the first part we describe our model, presenting the transport of heat by mixing. The formulation of the constrained MEP optimization problem is given (Sect. 2). Then, we compute the temperature, the specific energy, and the energy flux profiles. We give a physical interpretation of the effect of the constraint emerging from the positivity of the mass mixing coefficient. The impact of different expressions for energy is discussed (Sect. 3). A sensitivity test for the concentrations of O₃ and CO₂ is also performed. Finally, we discuss further works and objectives (Sect. 4). The computation of the geopotential is given in Appendix A and the resolution of the optimization problem is described in Appendix B.

2 Model

2.1 Vertical structure of the atmosphere

The atmosphere is divided into a column of N vertical layers. We work with prescribed pressure levels, so the elevation z depends on the temperature profile (see Appendix A). The CO₂, O₃, and water vapour profiles are fixed according to observations by McClatchey et al. (1972) and the ground is represented by a layer with a fixed surface albedo α. The atmosphere is supposed to be in hydrostatic equilibrium and is considered an ideal gas. The specific energy (energy per unit mass) in layer i – of mean elevation z_i, temperature T_i, and mixing ratio q_i (ratio between the mass of water vapour and total mass of the air for a given volume) – is the so-called moist static energy,

$$\begin{array}{}\text{(1)}& e_i=C_p{T}_i+g{z}_i+L{q}_i,\end{array}$$

where C_p=1005 J kg⁻¹ K⁻¹ is the heat capacity of the air, g=9.81 m s⁻² is the terrestrial acceleration due to gravity, and $L=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹ is the latent heat of vaporization.

We note ℛ_i, the net radiative energy input in layer i, by taking into account several effects: shortwave radiation, longwave radiation, reflection, and reabsorption. More explicitly,

$${\mathcal{R}}_i={\mathrm{SW}}_i+{\mathrm{LW}}_i={\mathrm{SW}}_{i\downarrow }-{\mathrm{SW}}_{i\uparrow }+{\mathrm{LW}}_{i\downarrow }-{\mathrm{LW}}_{i\uparrow },$$

where SW_i↓ is the downward radiative energy flux for short waves, SW_i↑ is the upward radiative energy flux for short waves, LW_i↓ is the downward radiative energy flux for long waves, and LW_i↑ is the upward radiative energy flux for long waves.

We use the code developed in Herbert et al. (2013) to compute the radiative budget. This model was developed to give a realistic description of the absorption properties of the more radiatively active constituents of the atmosphere while keeping a smooth dependence of the radiative flux with respect to the temperature profile. As suggested by the authors, this last requirement is important in the framework of a variational problem. The model is based on net exchange formalism (Dufresne et al., 2005), where the basic variables are the net exchange rates between each pair of layers instead of radiative fluxes.

In the longwave domain, the code decomposes the spectrum into 22 narrow bands and, in each band, it accounts for absorption by water vapour and carbon dioxide only. The absorption coefficient is computed using the statistical model of Goody (1952) with the data from Rodgers and Walshaw (1966). For the spatial integration, the diffusive approximation is performed with the standard diffusion factor $\mathit{\mu }=\mathrm{1}/\mathrm{1.66}$. Apart from the absorption data, given once and for all, the inputs of the model are the water vapour density, temperature profile, and carbon dioxide concentration. One may fix either absolute or relative humidity. In the shortwave domain, absorption by water vapour and ozone is accounted for by adapting the parameterization from Lacis and Hansen (1974). The input parameters for the model are the water vapour density and ozone density profiles, as well as surface albedo and solar constant. Clouds are not considered in the model. More details can be found in Herbert et al. (2013) and its supplementary material. The net radiative budget for the atmospheric layer i, ℛ_i, is given by summing over all terms involving the layer in question. In particular, ℛ_i is a function of all temperatures ${\left\{T_j\right\}}_{j=\mathrm{0},\mathrm{\dots },N}$ in the profile. Then,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\mathcal{R}}_i\left(T,q,{\mathrm{O}}_{\mathrm{3}},{\mathrm{CO}}_{\mathrm{2}},\mathit{\alpha }\right)=\\ {\displaystyle }& {\displaystyle }{\mathrm{SW}}_i\left(q,{\mathrm{O}}_{\mathrm{3}},\mathit{\alpha }\right)+{\mathrm{LW}}_i\left(T,q,{\mathrm{CO}}_{\mathrm{2}}\right).\end{array}\end{array}$$

In the previous equation, and in the following, T, q, O₃, and CO₂ will refer to complete profiles (i.e. $T=\mathit{\{}{T}_i{\mathit{\}}}_{i=\mathrm{0},\mathrm{\dots },N}$). Given that q (or $h=q/q_{\mathrm{s}}\left(T\right)$ fixed relative humidity), O₃, CO₂, and α are fixed in our model, we will only indicate the T dependence.

The vertical energy flux is represented by mixing between adjacent layers. Then, the net upward energy flux between layers i and i−1 writes as

$$\begin{array}{}\text{(3)}& F_i=m_i\left(e_{i-\mathrm{1}}-e_i\right),\end{array}$$

where m_i is a mixing coefficient that represents a mass per unit time and per unit surface. It is not (necessarily) homogeneous in all the column. We notice that m_i is typically the kind of coefficient that requires, at some point, a parameterization in usual climate models. Here m is obtained from the MEP procedure.

Taking into account the net radiative energy budget ℛ_i, the energy balance at the stationary state for the layer i reads as

$$\begin{array}{}\text{(4)}& F_i-F_{i+\mathrm{1}}+{\mathcal{R}}_i=\mathrm{0}.\end{array}$$

Figure 1 Discretization of an atmospheric column into N layers. The layer i – at temperature T_i, fixed pressure P_i, elevation z_i (which depends on the temperature profile), and mixing ratio q_i – has a specific energy $e_i=C_p{T}_i+g{z}_i+L{q}_i$; m_i is the mass mixing coefficient between layers i−1 and i which leads to the net energy flux $F_i=m_i(e_{i-\mathrm{1}}-e_i)$; ℛ_i is the net radiative energy budget in the layer i. The ground is represented by layer 0 with fixed surface albedo α.

2.2 Maximum entropy production with constraint

The principle of our model is to determine the fluxes F and temperatures T with the maximization of the entropy production. In the thermodynamics of diffusive processes, the entropy production is expressed as the sum of products of the fluxes with their associated thermodynamic forces. But we aim here at representing convection, not diffusion. We need to represent how the air parcels are mixed but only after (turbulent) convective motions. We can usefully separate this process into two steps. First, we assume a pseudo-adiabatic motion of an air parcel from one layer to another due to convection. This step is purely mechanical, without entropy production since we neglect viscous dissipation. The energy of the air parcel is conserved but its temperature and composition may change during the motion. If the air parcel was initially in layer i and goes to the layer i+1, the temperature of the air parcel becomes, by conservation of energy,

$$\begin{array}{}\text{(5)}& T_i^{\prime }=T_i+{\displaystyle \frac{\mathrm{1}}{C_p}}\left[g(z_i-z_{i+\mathrm{1}})+L(q_i-q_{i+\mathrm{1}})\right].\end{array}$$

Here, we have assumed that the water vapour concentration changes pseudo-adiabatically during the convection and not due to the mixing. This is not fully consistent since we do not impose water conservation. Secondly, the air parcel is mixed by diffusion with the ambient air in layer i+1 and transfers an amount of sensible heat per unit mass $C_p{T}_i^{\prime }$. At the same time, air parcels leave the layer i+1 for the layer i with a sensible heat per unit mass C_pT_i+1. So the net flux of sensible heat due to this process is

$$\begin{array}{}\text{(6)}& m_{i+\mathrm{1}}{C}_p(T_i^{\prime }-T_{i+\mathrm{1}})=m_{i+\mathrm{1}}(e_i-e_{i+\mathrm{1}})=F_{i+\mathrm{1}},\end{array}$$

where m_i+1 is the mass mixing coefficient between layers i and i+1. So the entropy production that results due to the sensible heat exchange between layers i and i+1 is the product of the flux F_i+1 with the thermodynamic force associated with the sensible heat only (i.e. the gradient of inverse temperature)

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{i+\mathrm{1}}=F_{i+\mathrm{1}}\left({\displaystyle \frac{\mathrm{1}}{T_{i+\mathrm{1}}}}-{\displaystyle \frac{\mathrm{1}}{T_i}}\right).\end{array}$$

By summing over all layers, and using the fact that $F_{N+\mathrm{1}}=\mathrm{0}$, we show that the total entropy production can be written as

$$\begin{array}{}\text{(8)}& \mathit{\sigma }=\sum _{i=\mathrm{0}}^N{\displaystyle \frac{\left(F_i-F_{i+\mathrm{1}}\right)}{T_i}}.\end{array}$$

In thermodynamics, more terms may contribute to entropy production such as pressure-volume work, or mixing. As for other MEPMs (Kleidon, 2010), we only retain the sensible heat exchange term. The geopotential and latent heat terms appear in the entropy production only as a result of our representation of convective transport, which is supposed to occur as a mechanically induced mass transport without entropy production. We can easily express the entropy production with temperatures σ(T) by using the energy balance in stationary state, $F_i-F_{i+\mathrm{1}}+{\mathcal{R}}_i\left(T\right)=\mathrm{0}$. Then the problem is usually solved in terms of temperature with a global constraint of energy conservation:

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

Given the form of the energy transport (Eq. 3), we here need to have additional constraints. Namely, the mass mixing coefficients m_i must be positive. It is then natural to solve the problem of flux by expressing the entropy production σ(F) and inequality constraints m_i≥0 with energy fluxes. Assuming that the relation ℛ(T) is invertible (Appendix B), we can formally write

$$\begin{array}{}\text{(10)}& F_{i+\mathrm{1}}-F_i={\mathcal{R}}_i\left(T\right)\iff T_i={\mathcal{R}}_i^{-\mathrm{1}}\left(F\right).\end{array}$$

This results in the following optimization problem with inequality constraints:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}\left.\underset{F_{\mathrm{1}},\mathrm{\dots },F_N}{max}\left(\sum _{i=\mathrm{0}}^N{\displaystyle \frac{F_i-F_{i+\mathrm{1}}}{{\mathcal{R}}_i^{-\mathrm{1}}\left(F\right)}}\right)\right|\exists m_i\ge \mathrm{0}\\ \\ \text{with}{F}_i=-m_i\left(e_i-e_{i-\mathrm{1}}\right)\end{array}\right\}.\end{array}$$

The constraint $F_i=-m_i\left(e_i-e_{i-\mathrm{1}}\right)$ with m_i≥0 naturally depends on the specific energy e used in the model. The later simply imposes the energy fluxes to be opposed to the energy gradient. We point out that the energy conservation is implicit in the flux formulation of the variational problem and does not need to be imposed as a constraint here.

3 Results

We have computed temperature, specific energy (energy per unit mass), and energy flux profiles for different prescribed atmospheric compositions from McClatchey et al. (1972) corresponding to the tropical, mid-latitude summer, mid-latitude winter, sub-arctic summer, and sub-arctic winter conditions. We work with a fixed relative humidity profile (ratio of the partial pressure of water vapour to the equilibrium vapour pressure of water at a given temperature). Typical values of surface albedo are used: α=0.1 for the tropical and mid-latitude conditions and α=0.6 for sub-arctic ones. The atmosphere is discretized in N=20 vertical levels.

3.1 The effect of the constraint

We investigate the effect of the following energy terms on the constraint:

• sensible heat C_pT;

• geopotential gz;

• latent heat for a water-vapour-saturated air Lq_s(T), where q_s is the mixing ratio at the saturation point. Since we work in pressure coordinates, it depends only on local temperature. This is a first attempt to take into account the effect of humidity without an explicit derivation of the humidity profile and water cycle. However, the radiative budget is still computed using a fixed standard relative humidity profile.

For illustration purposes, we can consider the case with only two layers. We note $F=m(e_{\mathrm{1}}-e_{\mathrm{2}})$ the net energy flux from layer 1 to layer 2, where m is the mass mixing coefficient between layers (cf. Fig. 2). In this simple case, the entropy production is written as $\mathit{\sigma }=F\left(\mathrm{1}/T_{\mathrm{2}}-\mathrm{1}/T_{\mathrm{1}}\right)$ and is limited by the constraint $m\ge \mathrm{0}⟺F(e_{\mathrm{1}}-e_{\mathrm{2}})\ge \mathrm{0}$. We interpret the different energy terms on this constraint as follow.

Figure 2 Energy exchanges between two layers of elevation z₁ and z₂ (z₂≥z₁), temperatures T₁ and T₂, and specific energy e₁ and e₂. We note $F=m(e_{\mathrm{1}}-e_{\mathrm{2}})$, the energy flux from layer 1 to layer 2, where m is the mass mixing coefficient between the two layers.

• e=C_pT: F≥0 if T₁≥T₂. The constraint simply imposes the energy transport from hot to cold regions.

• $e=C_{p}T+gz$: F≥0 if $T_{\mathrm{1}}\ge T_{\mathrm{2}}+g(z_{\mathrm{2}}-z_{\mathrm{1}})/C_p$. The geopotential gz limits the upward energy flux. We predict a warmer air at the bottom and a colder air at the top compared to the model with only sensible heat. Mechanically, we can see this as the expression of the fact that an air parcel from layer 1 may not have enough energy to move adiabatically to layer 2.

• $e=C_{p}T+gz+L{q}_{\mathrm{s}}$: F≥0 if $T_{\mathrm{1}}\ge T_{\mathrm{2}}+\left[g(z_{\mathrm{2}}-z_{\mathrm{1}})+L\left(q_{\mathrm{s}}\right(T_{\mathrm{2}})-q_{\mathrm{s}}(T_{\mathrm{1}}\left)\right)\right]/C_p$. Since q_s(T) is an increasing function, q_s(T₂)−q_s(T₁) has the same sign as T₂−T₁. The temperature gradient is usually negative (i.e. T₂≤T₁). Adding the latent heat at saturation makes the upward transport of energy less constrained. Consequently, the atmospheric temperature gradient weakens. Mechanically, we can see this as the fact that moist convection is easier than dry convection because of the transport of latent energy from the bottom to the top of the atmosphere.

Figure 3 Energy flux F, specific energy e, and temperature T for tropical atmospheric composition measured by McClatchey et al. (1972) and different expressions for energy on the constraint (e=C_pT, $e=C_{p}T+gz$, and $e=C_{p}T+gz+L{q}_{\mathrm{s}}$). The elevation is given in pressure level P. Results for the unconstrained model of Herbert et al. (2013) are represented. We also give the temperature profile corresponding to the measurements of McClatchey et al. (1972), labelled as “reference” for qualitative comparison.

3.2 General remarks

Various profiles are shown for tropical (Fig. 3) and sub-arctic winter (Fig. 4) conditions. The outputs of our constrained model are labelled by the energy terms taken into account in the constraint (e=C_pT, $e=C_{p}T+gz$, or $e=C_{p}T+gz+L{q}_{\mathrm{s}}$). The represented energy profiles are the energy corresponding to the constraint. For C_pT, we represent the profile e=C_pT, for C_pT+gz we represent $e=C_{p}T+gz$, and for $C_{p}T+gz+L{q}_s$ we represent $e=C_{p}T+gz+L{q}_s$. For e=C_pT, the specific energy is trivially more important for hot regions. For $e=C_{p}T+gz$, the geopotential adds energy to upper layers. For $e=C_{p}T+gz+L{q}_s$, the latent heat term adds energy to more humid layers. The results of the “unconstrained” model of Herbert et al. (2013) are represented with the associated thermal energy e=C_pT (rigorously speaking, the model is constrained by the global conservation of energy, but we will refer to it as unconstrained since no constraint is imposed on fluxes). Temperature profiles measured by McClatchey et al. (1972) are also represented for qualitative comparison and labelled as “reference”.

Figure 4 Energy flux F, specific energy e, and temperature T for sub-arctic winter atmospheric composition measured by McClatchey et al. (1972) and different expressions for energy on the constraint (e=C_pT, $e=C_{p}T+gz$, and $e=C_{p}T+gz+L{q}_{\mathrm{s}}$). The elevation is given in pressure level P. Results for the unconstrained model of Herbert et al. (2013) are represented. We also give the temperature profile corresponding to the measurements of McClatchey et al. (1972), labelled as “reference” for qualitative comparison.

For the unconstrained model, the energy flux is positive (i.e. upward) for the tropical (Fig. 3) and sub-arctic winter (Fig. 4) atmospheric compositions in all the column, despite the energy gradient inversion in the upper layers of the atmosphere. Therefore, the flux is in the same direction as the energy gradient in this region. This also corresponds to local negative entropy production and is not physically relevant. Consequently, the upward flux is overestimated, and the temperature gradient is weak. As discussed above (Sect. 2.1), the addition of the constraint m_i≥0 imposes the energy flux to be opposed to the specific energy gradient everywhere. As a consequence, the local entropy production is also constrained to be positive as in Ozawa and Ohmura (1997). But with an explicit account for the mass transport as explained above, we can account not only for sensible heat but more generally for moist static energy transfers in a convective column. If an energy flux in this direction is not favourable in matters of entropy production, it vanishes and we have stratification. When the geopotential term is considered (i.e. $e=C_{p}T+gz$ or $e=C_{p}T+gz+L{q}_{\mathrm{s}}$), we observe the following.

• An energy profile divided into three regions:

1. An unstable surface layer with a decreasing energy profile;

2. A neutral (slightly stable) mixed layer in the middle atmosphere with a vanishing energy gradient;

3. An inversion layer at the top of the atmosphere where the energy is increasing with elevation.

• A vanishing energy flux (stratification) in the upper part of the atmosphere (around P≃300 hPa for the tropics, Fig. 3, and P≃700 hPa for sub-arctic winter, Fig. 4).

We note that the thermal gradient is divided roughly by a factor of 2 when geopotential is considered, which gives a more realistic temperature profile.

3.3 Comparison between profiles

3.3.1 Model outputs for different climatic conditions

The constrained model is obviously sensitive to the water content. Considering $e=C_{p}T+gz$ or $e=C_{p}T+gz+L{q}_{\mathrm{s}}$ in the constraint gives approximately the same results for sub-arctic winter conditions (Fig. 4) (since q_s(T) is weak for low temperatures) while predictions differ for tropical conditions (Fig. 3) (where T and then q_s(T) are more important). We have verified that the influence of the surface albedo explains a large part of the temperature modification when we compare different climatic conditions.

3.3.2 Model output vs. reference profile

Before discussing the differences between the outputs and observations, we point out that this conceptual model does not take into account some important physical processes:

• Insolation is assumed to be constant, fixed at 1368∕4 W m⁻² for all conditions. But in reality, it varies with respect to seasons, diurnal cycle, and latitude due to the Earth's geometry and obliquity. So, the model does not take into account the variation in the radiative budget because of these geometrical factors;

• Horizontal energy fluxes are not considered in this 1-D description;

• The effect of clouds, which plays an important role in the absorption and emission of radiation (Dufresne and Bony, 2008), is not implemented in the radiative code.

Therefore, the aim of this study is not to give realistic values of temperature profiles nor vertical energy fluxes but to give a qualitative evaluation of the model. However, we can make some remarks on the comparison of our results to the reference temperature profiles. We observe that our model with $e=C_{p}T+gz+L{q}_{\mathrm{s}}$ provides better results for tropical conditions (Fig. 3), whereas the computed profiles are not so good for sub-arctic winter conditions (Fig. 4). Considering the previous remarks, one can explain the gap between our model and observations as follows. Constant insolation at the seasonal timescale is valid for the tropics, but it varies strongly for high latitudes. So our model is not adapted to represent a specific season at a high latitude like sub-arctic winter conditions.

Tropical regions are subject to strong vertical motion due to radiative heating. So horizontal energy fluxes are less important and the 1-D vertical description may be more adapted to this case. In contrast, radiative heating is less important for the Arctic (especially in winter), so convection is weaker. Then, the representation of horizontal energy fluxes is essential at high latitudes since they play a major role in heat transport from hot equatorial regions to cold poles. This probably explains why we underestimate temperature for high latitudes (Fig. 4).

3.4 Sensitivity to atmospheric composition

3.4.1 Ozone

When the influence of O₃ is not considered in the radiative budget, we observe small changes for specific energy and large changes for temperature in the stratosphere (cf. Fig. 5). Indeed, O₃ absorbs solar radiation at the top of the atmosphere, which induces heating of this region. It follows that the temperature in the high atmosphere is more important with ozone, and we even observe an inversion of the temperature gradient. It follows that less solar radiation heats the ground, resulting in smaller surface temperature. For our constrained model with $e=C_{p}T+gz+L{q}_{\mathrm{s}}$ including ozone, we observe a downward convective energy flux at the top of the atmosphere (Fig. 5). Ozone is therefore associated with heating from the top with an inversion of the temperature gradient and downward energy fluxes in the high atmosphere. This last effect only appears when both geopotential energy and O₃ are taken into account. When geopotential is not considered, the upward energy flux is so overestimated that the effect is undetectable (Figs. 3 and 4).

Figure 5 Energy flux F, specific energy e, and temperature T, with and without ozone, for the constrained model with $e=C_{p}T+gz+L{q}_{\mathrm{s}}$ and for the unconstrained model and for tropical atmospheric composition by McClatchey et al. (1972). Oc: with O₃, constrained. Nc: without O₃, constrained. Ou: with O₃, unconstrained. Nu: without O₃, unconstrained.

3.4.2 Carbon dioxide

We also have performed the classic experiment of doubling CO₂ concentration (Randall et al., 2007). The climate sensitivity, defined as the surface temperature differences between computations with [CO₂]=560 ppm and [CO₂]=280 ppm, is reported for the different atmospheric compositions in Table 1. Conventional models usually represent various processes like water vapour, ice-albedo, lapse rate, and cloud feedbacks. They play an important role in amplifying the climate sensitivity (Forster and Gregory, 2006). When comparing our values with the literature, we must keep in mind that our model does not represent all those feedbacks. The lapse rate feedback is taken into account. Water vapour feedback is partially represented in a crude way by fixing relative humidity (changes in temperature have an impact on water content and change the radiative budget) but there is no explicit representation of the hydrological cycle. It is technically possible to include the ice-albedo feedback in a MEPM (Herbert et al., 2011a), but this is not the case here. Clouds are not represented in the model so we will focus on the comparison between the constrained and unconstrained models using the same radiative scheme. Nevertheless, typical values of climate sensitivity for multi-model averages with only relevant feedbacks are given (Dufresne and Bony, 2008) for qualitative comparison. For fixed absolute water profile, our values are compared to Dufresne and Bony (2008) accounting only for the lapse rate feedback; for the fixed relative humidity, we compare them to Dufresne and Bony (2008) accounting for both lapse rate and water vapour feedbacks.

The sensitivity values computed here for the unconstrained model differ from Herbert et al. (2013). We have checked that it is only due to the fact that we use N=20 atmospheric layers here instead of N=9 in Herbert et al. (2013) (see Appendix B for the convergence of the algorithm with N). The climate sensitivity is higher for the constrained model than the unconstrained one (despite one exception with fixed absolute moisture for mid-latitude winter). This result may be interpreted as follow. If we start with a radiative forcing induced by a CO₂ doubling, it is the same for the two models since we fix identical atmospheric compositions and surface albedos. However, upward energy fluxes are limited for the constrained case which induces more important warming of the lower part of the atmosphere. The induced ground temperature increase is, therefore, more important for the constrained model. This effect is observable when we look at the perturbation of energy flux, specific energy, and temperature (Fig. 6). The constrained model provides more realistic values of sensitivity for fixed relative humidity, particularly for the tropics. Indeed, the sensitivity of 1.60 K computed in this case is closer to the literature value of 2.1±0.2 K (Dufresne and Bony, 2008) compared to the unconstrained model.

(Dufresne and Bony, 2008)(McClatchey et al., 1972)

Table 1 Climate sensitivity (warming, in K, of the surface due to a doubling of CO₂ concentration) of the constrained model with $e=C_{p}T+gz+L{q}_{\mathrm{s}}$, unconstrained model of Herbert et al. (2013), and literature (Dufresne and Bony, 2008). We give the values for different atmospheric compositions and for fixed absolute or relative water vapour profiles.

Figure 6 Differences in convective energy flux F, specific energy e, and temperature T between [CO₂]=560 ppm and [CO₂]=280 ppm for tropical atmospheric composition by McClatchey et al. (1972) represented for constrained model with e=C_pT, $e=C_{p}T+gz$, and $e=C_{p}T+gz+L{q}_{\mathrm{s}}$, and for the unconstrained model by Herbert et al. (2013).

4 Discussion

MEPMs are different from the usual GCMs or EMICs. Generally, atmospheric models are based on the following.

1. Kinematics: equations describing how the fluid moves.

2. Dynamics: equations describing why the fluid moves. They are based on Navier–Stokes equations linking the fluid acceleration to the forces.

3. Thermodynamics: energy budget equation involving dissipation, radiation, phases changes, etc.

4. An equation of state such as the ideal gas relation and approximations such as hydrostatic equilibrium to simplify the problem.

5. Closure hypothesis and/or parameterizations to represent sub-grid processes.

In usual climate models, energy transport is obtained after the computation of the velocity and other fields, whereas in MEPMs energy fluxes are computed implicitly without consideration of the dynamics. According to Dewar (2003) and Dewar (2009), MEP might be viewed as a statistical inference (such as the information theory interpretation of statistical physics by Jaynes, 1957) rather than a physical law. From this point of view, MEP allows drawing predictions (heat fluxes, specific energy, and temperature profiles) from the partial knowledge of the radiative budget and the energy content, while the effect of the sub-grid-scale processes is unknown. There are recent attempts to link MEP to other variational principles for dynamical systems or non-equilibrium statistical physics like the maximum of Kolmogorov–Sinaï entropy (Mihelich, 2015). Similar methods like the maximization of dynamical entropy or “Maximum Caliber” (Monthus, 2011; Dixit et al., 2018) may also be relevant to understand where MEP arises from a more formal point of view.

The gap between the model and observations is easily understood and explains why our 1-D description with constant insolation is more adapted for the tropics than arctic conditions. Further improvements are needed to solve these problems. Firstly, a more general 2- or 3-dimensional mass scheme transport is required. Secondly, it is formally possible to compute winds (Karkar and Paillard, 2015), moisture profiles, and humidity fluxes using MEP. Finally, we need to include a time dependence in our model for seasonal or diurnal cycles. A long-term objective might be to construct an SCM, with a limited number of adjustable parameters.

5 Conclusions

We have investigated the possibility of computing the vertical energy fluxes and temperature in the atmosphere using the MEP closure hypothesis into a simple climate model. The fluxes are then computed in an implicit way, which avoids tuning parameters. Contrary to some authors (Ozawa and Ohmura, 1997; Pujol and Fort, 2002; Herbert et al., 2013; Pascale et al., 2012), we have given a description of how the energy is transported. This paper provides the first attempt, to the best of the authors' knowledge, to introduce such a representation in a MEPM. Different energy terms can be considered: sensible heat, geopotential, and latent heat for a saturated profile. We have shown that this better energetic description of convection allows for obtaining more physically relevant temperature, specific energy, and energy flux profiles, still without any adjustable parameter for the dynamics. In particular, considering geopotential leads to stratification in the upper atmosphere and allows us to reproduce a temperature gradient closer to the observed one. We have investigated the sensitivity of the model when the atmospheric composition is modified. The results were compared to previous MEPMs and the literature. Our model is more sensitive to CO₂ than in Herbert et al. (2013) because the geopotential limits upward energy fluxes.

We hope that the present model may be helpful to construct SCMs with a reduced number of adjustable parameters.

Appendix A: Computing geopotential

We show here how to compute the geopotential and the dry static energy $e_i^{\mathrm{d}}=C_p{T}_i+g{z}_i$ as a function of temperature. One first writes

$$\begin{array}{}\text{(A1)}& z_i=z_i-z_{i-\frac{\mathrm{1}}{\mathrm{2}}}+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}\mathrm{\Delta }{z}_j,\end{array}$$

where $\mathrm{\Delta }{z}_j=z_{j+\frac{\mathrm{1}}{\mathrm{2}}}-z_{j-\frac{\mathrm{1}}{\mathrm{2}}}$ is the height of the layer j. So if ρ is the density of the air, R is the specific air constant and we assume the atmosphere is an ideal gas at hydrostatic equilibrium

$$\begin{array}{}\text{(A2)}& g\mathrm{\Delta }{z}_j=\underset{z_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{z_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}g\mathrm{d}z=-\underset{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}{\displaystyle \frac{\mathrm{d}p}{\mathit{\rho }}}=-\underset{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{\overset{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{\int }}RT{\displaystyle \frac{\mathrm{d}p}{p}}.\end{array}$$

Then, we can compute the mean elevation of a layer with two possible prescriptions:

• Isothermal layers (T=T_j in the integrand).

  $$\begin{array}{}\text{(A3)}& g\mathrm{\Delta }{z}_j=R{T}_j\ln \left({\displaystyle \frac{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}\right).\end{array}$$

  So the geopotential reads

  $$\begin{array}{}\text{(A4)}& g{z}_i=R\left[T_i\ln \left({\displaystyle \frac{p_{i-\frac{\mathrm{1}}{\mathrm{2}}}}{p_i}}\right)+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{T}_j\ln \left({\displaystyle \frac{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}\right)\right].\end{array}$$

• Dry isentropic layers ($T=T_j{\left(\frac{p}{p_j}\right)}^{\frac{R}{C_p}}$ in the integrand):

  $$\begin{array}{}\text{(A5)}& g\mathrm{\Delta }{z}_j=C_p{T}_j\left[{\left({\displaystyle \frac{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{p_j}}\right)}^{\frac{R}{C_p}}-{\left({\displaystyle \frac{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{p_j}}\right)}^{\frac{R}{C_p}}\right].\end{array}$$

  So the geopotential reads

  $$\begin{array}{}\text{(A6)}& g{z}_i=C_p\left[\begin{array}{l}{T}_i\left({\left({\displaystyle \frac{p_{i-\frac{\mathrm{1}}{\mathrm{2}}}}{p_i}}\right)}^{\frac{R}{C_p}}-\mathrm{1}\right)+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{T}_j\\ \left({\left({\displaystyle \frac{p_{j-\frac{\mathrm{1}}{\mathrm{2}}}}{p_j}}\right)}^{\frac{R}{C_p}}-{\left({\displaystyle \frac{p_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{p_j}}\right)}^{\frac{R}{C_p}}\right)\end{array}\right].\end{array}$$

In both cases, for imposed pressure levels, we obtain the following expression of the specific energy,

$$\begin{array}{}\text{(A7)}& e_i^{\mathrm{d}}\equiv \sum _{j=\mathrm{0}}^N\left(C_p{\mathit{\delta }}_{ij}+G_{ij}\right)T_j\equiv \sum _{j=\mathrm{0}}^N{E}_{ij}^{\mathrm{d}}{T}_j,\end{array}$$

where δ_ij is the Kronecker symbol and G_ij and $E_{ij}^{\mathrm{d}}$ are the coefficients of constant matrices G and E^d.

Appendix B: Resolution

In order to solve the optimization problem (Eq. 11), we express it in Lagrangian formalism, assuming strong duality holds (Boyd and Vandenberghe, 2004). We therefore search the critical points of the Lagrangian associated with this problem

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}\mathcal{L}& =\mathit{\sigma }-\sum _{i=\mathrm{1}}^N{\mathit{\mu }}_i{m}_i\\ & \text{with}{m}_i\ge \mathrm{0},{\mathit{\mu }}_i\ge \mathrm{0}\text{and}\\ & {\mathit{\mu }}_i{m}_i=\mathrm{0}\mathrm{for}i=\mathrm{1},\mathrm{\dots },N,\end{array}\end{array}$$

where ${\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\mu }}_N$ are Lagrange multipliers associated with the constraint (mass flux positivity). In order to formulate the problem in terms of energy fluxes F, we must express the inverse temperature $X=\mathrm{1}/T$ and the mass mixing m with F.

We use an iterative method to solve this non-linear optimization problem:

1. We linearize the radiative budget and specific energy around a given temperature profile.

2. The entropy production and constraints are then quadratic forms of energy fluxes that can be solved numerically.

3. We reiterate step 1 by linearizing around the temperature profile obtained in step 2 until convergence.

This is a rather standard procedure for optimization though there is no guarantee of finding the global solution in case of multiple local maxima. To overcome this issue, we start with various random initial temperature profiles that may lead to different local maxima. In the end, we retain only the best maximum, which is assumed to be the maximum of entropy production.

B1 Flux–temperature relation

The energy balance equation in stationary state can be written as follows:

$$\begin{array}{}\text{(B2)}& {\mathcal{R}}_i\left(X\right)+F_i-F_{i+\mathrm{1}}=\mathrm{0}.\end{array}$$

At each iteration, we linearize the radiative budget ℛ_i(X) around a reference temperature profile X⁰:

$$\begin{array}{}\text{(B3)}& {\mathcal{R}}_i\left(X\right)\simeq {\mathcal{R}}_i\left(X^{\mathrm{0}}\right)+\sum _{j=\mathrm{0}}^N{R}_{ij}(X_j-X_j^{\mathrm{0}}),\end{array}$$

where R is a square matrix of size N and ℛ(X⁰) is the radiative budget for the profile X⁰. If we assume R to be invertible, the energy flux can be computed with

$$\begin{array}{}\text{(B4)}& X_i=X_i^{\mathrm{0}}-\sum _{j=\mathrm{0}}^N{R}_{ij}^{-\mathrm{1}}\left(F_j-F_{j+\mathrm{1}}+{\mathcal{R}}_j\left(X^{\mathrm{0}}\right)\right).\end{array}$$

B2 Flux–mass relation

We first consider the dry static energy $e^{\mathrm{d}}=C_{p}T+gz$. Considering the atmosphere is an ideal gas at hydrostatic equilibrium, and for prescribed pressure levels, layer volume depends only on temperature. Therefore, the elevation of a layer is a function of temperatures of layers below only and we can express the energy of a layer as (see Appendix A)

$$\begin{array}{}\text{(B5)}& e_i^{\mathrm{d}}=\sum _{k=\mathrm{0}}^N{E}_{ik}^{\mathrm{d}}{T}_k,\end{array}$$

where E^d is a square, triangular matrix of size N. If we linearize around X⁰, one obtains

$$\begin{array}{}\text{(B6)}& e_i^{\mathrm{d}}\simeq \sum _{k=\mathrm{0}}^N{\displaystyle \frac{E_{ik}^{\mathrm{d}}}{X_k^{\mathrm{0}\mathrm{2}}}}(\mathrm{2}{X}_k^{\mathrm{0}}-X_k).\end{array}$$

Using Eq. (B4) gives the expression of energy as a function of flux

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{e}_i^{\mathrm{d}}\left(F\right)& \simeq \sum _{k=\mathrm{0}}^N{\displaystyle \frac{E_{ik}^{\mathrm{d}}}{X_k^{\mathrm{0}\mathrm{2}}}}\\ & \left(X_k^{\mathrm{0}}+\sum _{j=\mathrm{0}}^N{R}_{kj}^{-\mathrm{1}}\left(F_j-F_{j+\mathrm{1}}+{\mathcal{R}}_j\left(X^{\mathrm{0}}\right)\right)\right).\end{array}\end{array}$$

We also can take into account the latent heat for a water-vapour-saturated atmosphere. The mixing ratio at the saturation point, q_s, depends only on temperature T (in K) and pressure p (in Pa). It is given by the Bolton equation (Bolton, 1980):

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}& q_{\mathrm{s}}(T,p)={\displaystyle \frac{\mathrm{622.0}{h}_{\mathrm{s}}\left(T\right)}{p-h_{\mathrm{s}}\left(T\right)}}\\ & \mathrm{\text{with}}{h}_{\mathrm{s}}\left(T\right)=\mathrm{6.112}\exp \left({\displaystyle \frac{\mathrm{17.62}(T-\mathrm{273.15})}{T-\mathrm{30.03}}}\right),\end{array}\end{array}$$

where h_s is saturation vapour pressure of the mixture (in Pa). At fixed pressure, the moist static energy at saturation $e^{\mathrm{s}}=C_{p}T+gz+L{q}_{\mathrm{s}}$ of layers is only a function of temperatures. If we linearize around the profile X⁰,

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{e}_i^{\mathrm{s}}& =e_i^{\mathrm{d}}+L{q}_{\mathrm{s}}\left({\displaystyle \frac{\mathrm{1}}{X_i}}\right)\simeq e_i^{\mathrm{d}}+L{q}_{\mathrm{s}}\left({\displaystyle \frac{\mathrm{1}}{X_i^{\mathrm{0}}}}\right)\\ & -L{\left.{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial T}}\right|}_{\frac{\mathrm{1}}{X_i^{\mathrm{0}}}}{\displaystyle \frac{X_i-X_i^{\mathrm{0}}}{X_i^{\mathrm{0}\mathrm{2}}}}.\end{array}\end{array}$$

Then, we can use the same reasoning as for the dry static heat and replace the matrix E^d by E^s to consider the effect of latent energy for a saturated moisture profile. However, the radiative budget is still computed with reference water vapour profiles. In the following, e can represent e^d or e^s.

B3 Constraint

By multiplying both sides of

$$\begin{array}{}\text{(B10)}& F_i=-m_i\left(e_i-e_{i-\mathrm{1}}\right)\end{array}$$

by $\left(e_i-e_{i-\mathrm{1}}\right)$, we obtain

$$\begin{array}{}\text{(B11)}& F_i\left(e_i-e_{i-\mathrm{1}}\right)=-m_i{\left(e_i-e_{i-\mathrm{1}}\right)}^{\mathrm{2}}.\end{array}$$

So the constraint m_i≥0 is equivalent to

$$\begin{array}{}\text{(B12)}& {\mathit{\alpha }}_i\left(F\right)\equiv -F_i\left(e_i\left(F\right)-e_{i-\mathrm{1}}\left(F\right)\right)\ge \mathrm{0}.\end{array}$$

B4 Associated Lagrangian in flux space

Using the linearized energy budget (Eq. B4) and the constraint (Eq. B12), the problem (Eq. 11) is supposed to be equivalent to the search of critical points of the following Lagrangian

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}\mathcal{L}(F,\mathit{\mu })& =\mathit{\sigma }\left(F\right)-\sum _{i=\mathrm{1}}^N{\mathit{\mu }}_i{\mathit{\alpha }}_i\left(F\right)\\ & =\sum _{i=\mathrm{0}}^N{X}_i\left(F_i-F_{i+\mathrm{1}}\right)-\sum _{i=\mathrm{1}}^N{\mathit{\mu }}_i{\mathit{\alpha }}_i\end{array}\end{array}$$

$$\begin{array}{}\text{(B14)}& \begin{array}{rl}\simeq & \sum _{i=\mathrm{0}}^N\left(X_i^{\mathrm{0}}-\sum _{j=\mathrm{0}}^N{R}_{ij}^{-\mathrm{1}}\left(F_j-F_{j+\mathrm{1}}+{\mathcal{R}}_j\left(X^{\mathrm{0}}\right)\right)\right)\\ & \left(F_i-F_{i+\mathrm{1}}\right)-\sum _{i=\mathrm{1}}^N{\mathit{\mu }}_i{\mathit{\alpha }}_i\left(F\right),\end{array}\end{array}$$

while respecting the Karush–Kuhn–Tucker (KKT) conditions

$$\begin{array}{}\text{(B15)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathcal{L}}{\partial F_i}}& =\mathrm{0}\text{with}\left\{\begin{array}{l}{\mathit{\alpha }}_i\left(F\right)\ge \mathrm{0},\\ {\mathit{\mu }}_i\ge \mathrm{0},\end{array}\right.\\ & \text{and}{\mathit{\mu }}_i{\mathit{\alpha }}_i\left(F\right)=\mathrm{0}\mathrm{for}i=\mathrm{1},\mathrm{\dots },N.\end{array}\end{array}$$

The problem is solved numerically by using an interior point method (Boyd and Vandenberghe, 2004).

B5 Convergence of the algorithm

In practice, the algorithm may fail to find the global optimum for large N (≥50) but is robust for N≤40. As N increases, the algorithm converges rapidly to a solution (Fig. B1). The choice N=20 is a good compromise between computation time and resolution.

Figure B1 Energy flux F, specific energy e, and temperature T computed by our constrained model with $e=C_{p}T+gz+L{q}_{\mathrm{s}}$ for imposed tropical atmospheric composition measured by McClatchey et al. (1972), with various resolutions N.