The half-order energy balance equation – Part 1: The homogeneous HEBE and long memories

The original Budyko–Sellers type of 1D energy balance models (EBMs) consider the Earth system averaged over long times and apply the continuum mechanics heat equation. When these and the more phenomenological box models are extended to include time-varying anomalies, they have a key weakness: neither model explicitly nor realistically treats the conductive–radiative surface boundary condition that is necessary for a correct treatment of energy storage.

In this first of a two-part series, I apply standard Laplace and Fourier techniques to the continuum mechanics heat equation, solving it with the correct radiative–conductive boundary conditions and obtaining an equation directly for the surface temperature anomalies in terms of the anomalous forcing. Although classical, this equation is half-ordered and not integer-ordered: the half-order energy balance equation (HEBE). A quite general consequence is that although Newton's law of cooling holds, the heat flux across surfaces is proportional to a half-ordered (not first-ordered) time derivative of the surface temperature. This implies that the surface heat flux has a long memory, that it depends on the entire previous history of the forcing, and that the temperature–heat flux relationship is no longer instantaneous.

I then consider the case in which the Earth is periodically forced. The classical case is diurnal heat forcing; I extend this to annual conductive–radiative forcing and show that the surface thermal impedance is a complex valued quantity equal to the (complex) climate sensitivity. Using a simple semi-empirical model of the forcing, I show how the HEBE can account for the phase lag between the summer maximum forcing and maximum surface temperature Earth response.

In Part 2, I extend all these results to spatially inhomogeneous forcing and to the full horizontally inhomogeneous problem with spatially varying specific heats, diffusivities, advection velocities, and climate sensitivities. I consider the consequences for macroweather (monthly, seasonal, interannual) forecasting and climate projections.

1 Introduction

Ever since Budyko (1969) and Sellers (1969) proposed a simple model describing the exchange of energy between the Earth and outer space, energy balance models (EBMs) have provided a straightforward way of understanding past, present, and possible future climates. The models usually have either zero or one spatial dimension, respectively representing the globally or latitudinally averaged meridional temperature distribution (for a review, see McGuffie and Henderson-Sellers, 2005, and North and Kim, 2017).

The fundamental EBM challenge is to model the way that imbalances in incoming shortwave and outgoing longwave radiation are transformed into changes in surface temperatures. In an energy-balanced climate state, the vertical flux imbalances are transported horizontally. Here I am primarily interested in the anomalies with respect to this state. When an external flux (forcing) is added, some of this anomalous imbalance is radiated to outer space, while some is converted into sensible heat and conducted into (or out of) the subsurface. This latter flux accounts for both energy storage and surface temperature changes as well as attendant changes in longwave emissions. EBMs avoid explicit treatment of this critical surface boundary condition, treating it phenomenologically in ways that are flawed; in this two-part paper, I show how they can easily be improved with significant benefits: first, the (idealized) homogeneous case (Part 1) and then the general horizontally inhomogeneous (2D) case (Part 2; Lovejoy, 2021).

First, consider box EBMs with zero horizontal dimensions as a model of the mean Earth temperature. These are based on two distinct assumptions, namely (a) that the rate that heat (S) is exchanged between the Earth and outer space (dS $/$ dt) is proportional to the difference between the surface temperature (T) and its long-term equilibrium value (T_eq), dS $/$ dt $\propto (T_{\mathrm{eq}}-T)$ (Newton's law of cooling, NLC), and (b) that this rate is also proportional to the rate of change of surface temperature: dS $/$ dt∝ dT $/$ dt. Budyko–Sellers models are on firmer ground: they start with the basic continuum mechanics heat equation with advective and diffusive heat transport. Yet they have no vertical coordinate and are therefore unable to correctly treat the surface conduction–radiation–energy storage issue. By restricting explicit treatment of energy transport to the horizontal, they resort to the ad hoc assumption that the vertical flux imbalances are redirected horizontally and meridionally. The original Budyko–Sellers models were of time-independent climate states, and there was no energy storage at all: the radiative imbalances were completely redirected. While this approximation may be reasonable for these long-term states, they became problematic as soon the original models were extended to include temporal variations (Dwyers and Petersen, 1975). While these time-varying extensions implicitly allow for subsurface energy storage, this implicit treatment is both unnecessary and unsatisfactory.

The basic physical problem is that anomalous radiative flux imbalances partly lead to heat conduction fluxes into the subsurface and partly to changes in longwave radiative fluxes. The part conducted into the subsurface is stored and may re-emerge, possibly much later. Starting with the heat equation, realistic and mathematically correct treatments involve the introduction of a vertical coordinate and the use of conductive–radiative surface boundary conditions (BCs). If one considers the horizontally homogeneous 3D problem in a semi-infinite medium with these mixed BCs and linearized longwave emissions, the problem is classical and can be straightforwardly solved using Laplace and Fourier techniques. Mathematically it turns out that the key is the surface layer that defines the surface vertical temperature gradient. The influence of the subsurface is only over a thin layer of the order of a few diffusion depths where most of the energy storage occurs. This depth depends on the specific heat per volume and the diffusivity, and it is estimated to be typically of the order of 100 m for the ocean (depending on its turbulent diffusivity) and less over land (see Appendix A, Part 2).

The exact treatment of this homogeneous problem confirms that Newton's law of cooling holds but shows that the classical box model relation between heat flux and the surface temperature is wrong: symbolically the correct relation is $\mathrm{d}S/\mathrm{d}t\propto d^{h}T/\mathrm{d}{t}^h$ with $h=\mathrm{1}/\mathrm{2}$ – not the phenomenological value h= 1. Physically, these fractional derivatives are simply convolutions, and in the Fourier domain they are power-law filters, in this case involving power-law storage (hence “memories”). The corresponding half-order energy balance equation (HEBE) has qualitatively much stronger storage than the short exponential memories associated with the standard integer-ordered (h= 1) box model derivatives.

Half-order derivatives have appeared in heat and diffusion problems since at least the 1960s (Meyer, 1960; Oldham and Spanier, 1972; Oldham, 1973; Oldham and Spanier, 1974). An equation mathematically identical to the homogeneous h= $\mathrm{1}/\mathrm{2}$ special case of the fractional energy balance equation (FEBE) was derived by Oldham (1973) as a short-time approximation of electrolyte diffusion in a spherical geometry, and Oldham and Spanier (1974) anticipated the present application by noting that half-order derivatives can be applied to “not one but an entire class of boundary value problems”. Later, half-order derivatives were developed by Babenko (1986) and have been regularly exploited in engineering heat transfer problems (see, e.g., Sierociuk et al., 2013, 2015, and references therein). The method is probably not more generally known since most applications are with fairly standard heat flux boundary conditions and other more familiar techniques can also be used.

More generally, fractional derivatives and their equations (Podlubny, 1999) have a history going back to Leibniz in the 17th century, and their development has exploded in the last decades (for books on the subject, see, e.g., Miller and Ross, 1993; Podlubny, 1999; Hilfer, 2000, West et al., 2003; Tarasov, 2010; Klafter et al., 2012; Klafter et al., 2012, Baleanu et al., 2012; Atanackovic et al., 2014).

Interestingly, the explicit or implicit application of fractional derivatives to model the Earth's temperature – and more recently the energy budget – has several antecedents arising from the wide range of spatial scaling symmetries of atmospheric fields respected by the fluid equations, models, and (empirically) the atmospheric fields themselves (see the reviews in Lovejoy and Schertzer, 2013; Lovejoy, 2019a). Since this includes the velocity field – whose spatial scaling implies scaling in time – it implies that power laws should be more realistic than exponentials. At first, this led to power-law climate response functions (CRFs) (Rypdal, 2012; van Hateren, 2013; Rypdal and Rypdal, 2014; Rypdal et al., 2015; Hebert, 2017; Hébert et al., 2021). However, without truncations, pure power-law CRFs lead to divergences: the “runaway Green's function effect” (Hébert and Lovejoy, 2015), whereby a model is unstable to infinitesimal step function increases in forcing and the equilibrium climate sensitivity is infinite. These can be tamed by a high-frequency truncation (Hebert, 2017; Hébert et al., 2021) or avoided by constraining forcings to return to zero (Rypdal, 2016; Myrvoll-Nilsen et al., 2020).

However, Lovejoy (2019a) and Lovejoy et al. (2021) argued that it is not the CRF itself but rather the Earth's heat storage mechanisms that respect the scaling symmetry. This hypothesis implies that the corresponding storage (the derivative term) in the energy balance equation (EBE) is of fractional rather than integer order: the fractional energy balance equation (FEBE). Denoting the order of the derivative term in the equation by h, it was empirically shown that if the derivative was of order h≈ 0.4–0.5 (rather than the classical EBE value h= 1), it could account for both the low-frequency multidecadal memory (Hebert, 2017; Hébert et al., 2021) needed for climate projections and the high-frequency macroweather (i.e., the regime at timescales longer than the lifetime of planetary structures, here monthly to decadal) memory needed for monthly, seasonal, and annual macroweather forecasts (Lovejoy, 2015; Lovejoy et al., 2015; Del Rio Amador and Lovejoy, 2019, 2021a, b). Indeed, the FEBE CRF can be used directly to make climate projections that are compatible with the Coupled Model Intercomparison Project 5 (CMIP5) multi-model ensemble mean projections but with substantially smaller uncertainties (Procyk et al., 2020). Finally, it is possible to generalize the classical (3D) continuum equation to the fractional heat equation from which the (inhomogeneous, 2D) FEBE governs the surface temperature (work in progress).

In spite of empirical and theoretical support, the FEBE is essentially a phenomenological global model; in this paper I show how – at least for the h= $\mathrm{1}/\mathrm{2}$ special case – it can be placed on a firmer theoretical basis while simultaneously extending it to two spatial dimensions. The model is for macroweather temperature anomalies, i.e., at timescales longer than the lifetimes of planetary structures, typically 10 d. Following Budyko and Sellers, the system averaged over weather scales is considered to be a continuum, justifying the application of the continuum mechanics heat equation. The starting point is thus the same as the classical EBMs: radiative, advective, and conductive heat transport using the standard continuum mechanics energy equation. Also, following the classical approaches, the longwave black-body radiation is treated in its linearized form.

This work is divided into two parts. The first part is classical; it focuses on the homogeneous heat equation, pointing out the consequence that with semi-infinite geometry (depth) and with (realistic) conductive–radiative boundary conditions, the surface temperature satisfies the homogeneous HEBE. I relate this to the usual box models, Budyko–Sellers models, and classical diurnal heating models including the notions of thermal admittance and impedance as well as complex climate sensitivities that are useful in understanding the annual cycle. I underscore the generality of the basic (long-memory) storage mechanism. The second part extends this work to the horizontal, first to the homogeneous case (but with inhomogeneous forcing, including a direct comparison with the classical latitudinally varying 1D Budyko–Sellers model on the sphere) and then – using Babenko's method – to the general inhomogeneous case. Part 2 also contains several appendices that discuss empirical parameter estimates, spatial statistics useful for empirical orthogonal functions, and understanding the horizontal scaling properties as well as the changes needed to account for spherical geometry.

2 The transport equations

2.1 Conductive and advective heat fluxes

In most of what follows, the Earth's spherical geometry plays no role, and I use Cartesian coordinates with the z axis pointing upwards and horizontal coordinates x= (x, y) (however, in Sect. 2.3 in Part 2 and Appendix C of Part 2, I treat the latitudinally varying case on a sphere). The horizontal is essentially the same as in the Budyko–Sellers model: horizontal diffusive and advective heat fluxes are atmospheric column averages lying on the surface (z= 0). What is new is the treatment of the vertical with radiative and conductive fluxes crossing the surface either into the subsurface (downward, the negative z direction where it can propagate to −∞) or to outer space (upward, z > 0) so that heat is effectively stored in the half-volume ($x,y,z\mathit{<}\mathrm{0}$). Although in principle this means that the entire semi-infinite region z≤ 0 is modeled, ultimately only the vertical surface temperature derivative is needed and this is well defined as long as the surface layer is of the order of a few diffusion depths (tens or hundreds of meters). Later, I show that the main equations only explicitly depend on the local relaxation times and climate sensitivities; the vertical and horizontal transport details are only implicit. Finally, the fields are assumed to be in the macroweather regime; i.e., they have been averaged over the weather–macroweather transition scale (about 10 d) or longer and possibly for tens or hundreds of kilometers in space (the space–time limits are not yet clear). Since 10 d is the typical lifetime of planetary atmospheric structures, much of the actual turbulent atmospheric transport processes are averaged out, giving some justification to the parametrization. Future work is needed to clarify several foundational issues.

Let us start with energy transport by diffusion: Fick's law ${\mathit{Q}}_{\mathrm{d}}=-\mathit{\rho }c\mathit{\kappa }\mathrm{\nabla }T$, where Q_d is the diffusive heat flux vector, κ is the thermal diffusivity, ρ the density, c the specific heat, and T($\mathit{x},z,t$) the temperature. Following standard energy balance models, use eddy diffusivities that are different in the horizontal (h) and vertical (v) as κ_h(x) and κ_v(x):

$$\begin{array}{}\text{(1)}& Q_{\mathrm{d}}=-\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}T-\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T}{\partial z}}\widehat{z};{\mathrm{\nabla }}_{\mathrm{h}}={\displaystyle \frac{\partial }{\partial x}}\widehat{x}+{\displaystyle \frac{\partial }{\partial y}}\widehat{y}.\end{array}$$

The circumflex indicates unit vectors. To include advection, consider the heat equation for a fluid in a horizontal velocity field v_h:

$$\begin{array}{}\text{(2)}& c\mathit{\rho }{\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}=-\mathrm{\nabla }\cdot {\mathit{Q}}_{\mathrm{d}};{\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}={\displaystyle \frac{\partial T}{\partial t}}+{\mathit{v}}_{\mathrm{h}}\cdot \mathrm{\nabla }T,\end{array}$$

where D $/$ Dt is the advective derivative. The heat equation is therefore

$$\begin{array}{}\text{(3)}& c\mathit{\rho }{\displaystyle \frac{\partial T}{\partial t}}=-c\mathit{\rho }{\mathit{v}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{\mathrm{h}}T+{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}T\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T}{\partial z}}\right).\end{array}$$

If the volumetric specific heat (cρ) is constant, using the continuity equation $\mathrm{\nabla }\cdot \left(c\mathit{\rho }{\mathit{v}}_{\mathrm{h}}\right)=\mathrm{0}$, one may write

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& c\mathit{\rho }{\displaystyle \frac{\partial T}{\partial t}}=-\mathrm{\nabla }\cdot ({\mathit{Q}}_{\mathrm{a}}+{\mathit{Q}}_{\mathrm{d}});\\ & {\mathit{Q}}_{\mathrm{a}}=c\mathit{\rho }{\mathit{v}}_{\mathrm{h}}(T-T_{\mathrm{0}});{\mathit{Q}}_{\mathrm{d}}=-\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}T-\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T}{\partial z}}\widehat{z}.\end{array}\end{array}$$

Q_a is the advective heat flux and T₀ is a constant reference temperature (it disappears when the divergence is taken). This is the classical fluid heat equation; it can readily be verified that it conserves energy (integrate both sides over a volume and then use the divergence theorem). κ_h (x), κ_v(x), and v_h(x) are taken to be independent of t and z; they are part of the climate state and are empirically determined so as to reproduce the time-independent climate temperature distribution. In future work, they could be given their own time-varying anomalies.

2.2 Radiative heat fluxes

At the surface, there is an incoming energy flux R_↓,

$$\begin{array}{}\text{(5)}& R_{\downarrow }(\mathit{x},t)=Q_{\mathrm{0}}\left(\mathit{x}\right)+F(\mathit{x},t),\end{array}$$

where F is the anomalous forcing and Q₀(x) is the local solar radiation:

$$\begin{array}{}\text{(6)}& Q_{\mathrm{0}}\left(\mathit{x}\right)=QS\left(\mathit{x}\right)a\left(\mathit{x}\right).\end{array}$$

Q is the mean top-of-the-atmosphere flux (≈ 341 W m²), S(x) is the dimensionless local solar constant with local co-albedo a(x) (in the notation of North and Kim, 2017), and the time-dependent part of the radiative balance is specified by the additional incoming energy flux, the “forcing” F(x, t). Although in this paper I mostly ignore temporal albedo variations (see, however, Sect. 3.3), they are important for studying temperature–albedo feedbacks and climate transitions. If needed, even if they include a (potentially nonlinear) temperature dependence, they are easy to incorporate. For example, they could be included in F by using $a\left(\mathit{x},t\right)=a_{\mathrm{0}}\left(\mathit{x}\right)+a_{\mathrm{1}}\left(\mathit{x},t,T\left(\mathit{x},t\right)\right)$ in place of a(x) in Eq. (6) and $F\left(\mathit{x},t\right)=F_{\mathrm{0}}\left(\mathit{x},t\right)+QS\left(\mathit{x}\right)a_{\mathrm{1}}\left(\mathit{x},t,T\left(\mathit{x},t\right)\right)$ in place of F in Eq. (5).

As usual, F(x, t) includes solar, volcanic, and anthropogenic forcings. However, since macroweather includes random internal variability, F(x, t) also includes a stochastic internal variability component. Finally, for macroweather scales shorter than a year, F could also include the annual cycle and therefore possible cyclical albedo variations due to seasonally varying cloudiness (Sect. 3.3). Alternatively, T and F can be deseasonalized in the usual way to yield standard monthly climate “normals” so that the mean anomalies are zero over the climate normal reference period.

$R_{\downarrow }\left(\mathit{x},t\right)$ is partially balanced by the outgoing $R_{\uparrow }\left(\mathit{x},t\right)$ that depends on the surface temperature and the effective emissivity ε(x):

$$\begin{array}{}\text{(7)}& R_{\uparrow }\left(\mathit{x},t\right)=\mathit{\sigma }\mathit{\epsilon }\left(\mathit{x}\right)T{\left(\mathit{x},\mathrm{0},t\right)}^{\mathrm{4}},\end{array}$$

where σ is the Stefan–Boltzmann constant. The R_↓ and R_↑ imbalance drives the system, and it implies that heat diffuses across the surface, which is the top boundary condition needed to solve Eq. 3 for T(x, z < 0, t):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\left.\left(\mathit{\sigma }\mathit{\epsilon }\left(\mathit{x}\right)T(\mathit{x},z,t{)}^{\mathrm{4}}+\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}(\mathit{x}){\displaystyle \frac{\partial T(\mathit{x},z,t)}{\partial z}}\right)\right|}_{z=\mathrm{0}}\\ & =Q_{\mathrm{0}}\left(\mathit{x}\right)+F(\mathit{x},t).\end{array}\end{array}$$

The derivative term ${\left.\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}\partial T/\partial z\right|}_{z=\mathrm{0}}=Q_{\mathrm{s}}$ is the conductive (sensible) heat flux across the surface into the Earth; see Fig. 1. The radiative fluxes thus impose a “mixed” conductive–radiative boundary condition involving both T and $\partial T/\partial z$ (they are a special case of “Robin” boundary conditions; Hahn and Ozisk, 2012). If we add the initial condition $T(\mathit{x},z,t=\mathrm{0})=\mathrm{0}$ (or later, $T(\mathit{x},z,t=-\mathrm{\infty })=\mathrm{0}$ and the Dirichlet boundary condition at great depth $T\left(\mathit{x},z=-\mathrm{\infty },t\right)=\mathrm{0}$ and assume that the system is periodic or infinite in the horizontal, then, in principle, these are enough to determine the temperature for T(x, z < 0, t > 0) (or eventually, $T(\mathit{x},z,t=-\mathrm{\infty })=\mathrm{0}$). Instead of avoiding this conductive–radiative BC below I show how it directly yields an equation for the surface temperature.

Figure 1 A schematic diagram showing the correct 3D energy balance equations with conductive–radiative surface boundary conditions. Q_s is the heat flux across the surface into the subsurface, and S is the energy stored in the subsurface per unit surface area. The picture illustrates the thin surface layer (whose thickness is of the order of the diffusion depth l_v with relaxation time τ; Eq. 20) in which the radiative exchanges between the Earth and outer space occur.

2.3 The climatological and anomaly fields

Now decompose the heat flux and temperature into time-independent (climatological) and time-varying (anomaly) components: Q_c, T_c, Q, and T. As usual, linearize the outgoing black-body radiation, although do so around the spatially varying surface temperature T_c(x,z= 0) (i.e., not the global average temperature), which yields spatially varying coefficients:

$$\begin{array}{}\text{(9)}& R_{\uparrow }\left(T_{\mathrm{c}}\left(\mathit{x},\mathrm{0}\right)+T\left(\mathit{x},\mathrm{0},t\right)\right)\approx R_{\uparrow }\left(T_{\mathrm{c}}\left(\mathit{x},\mathrm{0}\right)\right)+{\displaystyle \frac{T\left(\mathit{x},\mathrm{0},t\right)}{s\left(\mathit{x}\right)}},\end{array}$$

where T_c+T is the actual temperature. With climate sensitivity, it yields

$$\begin{array}{}\text{(10)}& s\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\sigma }\mathit{\epsilon }\left(\mathit{x}\right)T_{\mathrm{c}}{\left(\mathit{x},\mathrm{0}\right)}^{\mathrm{3}}}}.\end{array}$$

Since typical macroweather temperature anomalies are only a few degrees, the black-body emission is quite linear with the temperature anomaly. However, due to feedbacks, the proportionality coefficient – the climate sensitivity – as estimated in Eq. (10) is not accurate; below, simply consider s(x) to be an empirically determined function of position.

The incoming radiation at the location x drives the system. The radiative imbalance ΔR going into the subsurface is therefore equal to the conductive flux Q_s into the surface; it specifies the conductive–radiative surface boundary condition for T_c and the anomalies T:

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }R=Q_{\mathrm{s}};\mathrm{\Delta }R=R_{\downarrow }-R_{\uparrow };Q_{\mathrm{s}}=-Q_{\mathrm{d},z},\end{array}$$

where Q_d,z is the (upward) vertical component of the heat flux at the surface given by Fick's law: $Q_{\mathrm{d},z}=-\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\left.\frac{\partial T}{\partial z}\right|}_{z=\mathrm{0}}$. The conductive–radiative surface boundary conditions for the time-independent climate and anomaly temperatures is therefore

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \left(R_{\uparrow }\left(T_{\mathrm{c}}\right(\mathit{x},z\left)\right)+\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{c}}(\mathit{x},z)}{\partial z}}\right){|}_{z=\mathrm{0}}=Q_{\mathrm{0}}\left(\mathit{x}\right)\\ & \left({\displaystyle \frac{T(\mathit{x},z,t)}{s}}+\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T(\mathit{x},z,t)}{\partial z}}\right){|}_{z=\mathrm{0}}=F(\mathit{x},t),\end{array}\end{array}$$

where s, ρ, c, and κ are all presumed to be functions of x. Note that the conductive heat flux is a sensible heat flux; the boundary condition involves its vertical component that represents heat stored in the subsurface. While Eqs. (11) and (12) involve the vertical temperature derivative at the surface (i.e., over an infinitesimal layer), l_v= sρcκ_v defines the diffusion depth (typically ≈ 10–100 m in thickness; see Part 2) so that physically the model need only be realistic over this fairly shallow depth where most of the (anomalous) heat is stored.

Now, in the temperature equation (Eq. 3), replace T by T_c+T. The equation for the time-independent climate part is

$$\begin{array}{}\text{(13)}& \begin{array}{rl}c\mathit{\rho }{\displaystyle \frac{\partial T_{\mathrm{c}}}{\partial t}}& =\mathrm{0}=-\mathit{\rho }c{\mathit{v}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{\mathrm{h}}{T}_{\mathrm{c}}+{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}{T}_{\mathrm{c}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{c}}}{\partial z}}\right).\end{array}\end{array}$$

And the equation for the time-varying anomalies is

$$\begin{array}{}\text{(14)}& c\mathit{\rho }{\displaystyle \frac{\partial T}{\partial t}}=-\mathit{\rho }c{\mathit{v}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{\mathrm{h}}T+{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}T\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T}{\partial z}}\right).\end{array}$$

These equations must now be solved using boundary conditions in Eqs. (11) and (12) for T_c, T, and $T_{\mathrm{c}}=T=$ 0 at $z=-\mathrm{\infty }$ (all t), as well as $T(\mathit{x},z,t=\mathrm{0})$ =0 (or see below, $T(\mathit{x},z,t=-\mathrm{\infty })=\mathrm{0})$.

The separation into one equation for the time-invariant climate state and another for the time-varying anomalies is done for convenience. As long as the outgoing longwave radiation is approximately linear over the whole range of temperatures (as is commonly assumed in EBMs), this division involves no anomaly smallness assumptions or assumptions concerning their time averages; the choice of the reference climate depends on the application. Below, I choose anomalies defined in the standard way (although not necessarily with the annual cycle removed; Sect. 3.3); this is adequate for monthly and seasonal forecasts as well as 21st-century climate projections. However, a different choice might be more appropriate for modeling transitions between different climates including possible chaotic behaviors.

2.4 The climatological temperature distribution and Budyko–Sellers models

In order to simplify the problem, starting with Budyko (1969) and Sellers (1969), the usual approach to obtaining T_c is somewhat different. First, the climatological temperature field is only defined at z= 0, i.e., T_c(x) =T_c(x, 0). Without a vertical coordinate, the climatological radiative imbalance $Q_{\mathrm{0}}\left(\mathit{x}\right)-R_{\uparrow }\left(T_{\mathrm{c}}\right(\mathit{x}\left)\right)$ no longer forces the system via the vertical surface derivative (Eq. 11); instead, the imbalance is conventionally redirected in the meridional direction away from the Equator (Fig. 2).

Figure 2 A schematic diagram showing the Budyko–Sellers 1D energy balance equation obtained by latitudinal averaging and by redirecting the vertical imbalance away from the Equator.

To see how this works, return to Eq. (4) for the climatological component and put $\frac{\partial }{\partial z}=\mathrm{0}$:

$$\begin{array}{}\text{(15)}& {\mathit{Q}}_{\mathrm{c}}\left(\mathit{x}\right)={\mathit{Q}}_{\mathrm{c},\mathrm{a}}\left(\mathit{x}\right)+{\mathit{Q}}_{\mathrm{c},\mathrm{d}}\left(\mathit{x}\right)+\mathrm{sign}\left(y\right)\left(Q_{\mathrm{0}}\right(\mathit{x})-R_{\uparrow }(T_{\mathrm{c}}\left(\mathit{x}\right)\left)\right)\widehat{y}.\end{array}$$

In this formulation, one usually uses the latitude angle instead of the meridional coordinate y (see Part 2, Sect. 2.3). The direction of the redirected vertical flux is always away from the Equator (y= 0, and hence sign(y)); in any event, zonal fluxes will cancel when averaged over latitudinal bands.

The usual Budyko–Sellers type of models then average Q_c over lines of constant latitude, yielding a 1D model:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{\mathit{Q}}}_{\mathrm{c}}\left(y\right)& =\left(\mathit{\rho }c\right(v_y{\stackrel{\mathrm{‾}}{T}}_{\mathrm{c}}-{\mathit{\kappa }}_{\mathrm{h}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{T}}_{\mathrm{c}}}{\partial y}})\\ & +\mathrm{sign}\left(y\right)\left(Q_{\mathrm{0}}\right(y)-R_{\uparrow }({\stackrel{\mathrm{‾}}{T}}_{\mathrm{c}}\left)\right))\widehat{y},\end{array}\end{array}$$

where the overbar indicates averaging over all longitudes, x.

In the more popular Seller's version, the basic horizontal transport is due to the eddy thermal diffusivity, the κ_h term. There may also be a small advection velocity v; it is not considered to be a true physical velocity but only an ad hoc parameter needed to prevent κ_h from being negative (Sellers, 1969). The standard presentation (North et al., 1981; North and Kim, 2017) avoids the problem by using the diffusivity as in Sect. 3.1. The horizontal eddy diffusivity κ_h is often taken as the sum of contributions from water, water vapor, and air. In the pure Budyko version, there is no eddy diffusivity, and the heat flux is assumed to be proportional to the temperature difference with respect to a reference (e.g., mean) value $\left(\mathit{Q}{)}_y\propto (T-T_{\mathrm{0}})\right)$. Comparing this with Eq. (4) for Q_a, we see that this implies that Budyko horizontal heat fluxes are purely advective.

The final step to obtaining the energy equation is to take the divergence:

$$\begin{array}{}\text{(17)}& \mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{Q}}}_{\mathrm{c}}={\displaystyle \frac{\partial ({\stackrel{\mathrm{‾}}{\mathit{Q}}}_{\mathrm{c}}{)}_y}{\partial y}}=-\mathit{\rho }c{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{T}}_{\mathrm{c}}}{\partial t}}.\end{array}$$

Budyko and Sellers only considered the time-independent case and obtained

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left({\stackrel{\mathrm{‾}}{\mathit{Q}}}_{\mathrm{c}}\right(y){)}_y}{\partial y}}=\mathrm{0}\\ & {\stackrel{\mathrm{‾}}{\mathit{Q}}}_{\mathrm{c}}=\mathrm{const}.\end{array}\end{array}$$

By appropriately choosing a reference temperature (usually the global average), the constant can be adjusted for convenience. Somewhat later, Dwyers and Petersen (1975) considered the time-independent case (Eq. 17), which is second order in y. Subsequently, the model has been widely used for studying different past and future climates as well as the corresponding transitions. Note that the $\mathit{\rho }c\frac{\partial T_{\mathrm{c}}}{\partial t}$ term corresponds to energy storage; in the time-independent case there is no storage.

3 The classical origin of the fractional operators: conductive–radiative boundary conditions in a semi-infinite domain

3.1 The zero-dimensional homogeneous heat equation

3.1.1 The key parameters

No matter how the climate temperature equation is solved, the equation for the time-dependent anomaly temperature remains as in Eq. (14). I now rewrite it in a way that brings out the critical mathematical properties. Since ρc and κ_v are taken to be only functions of x, Eq. (14) can be rewritten as

$$\begin{array}{}\text{(19)}& \left({\displaystyle \frac{\partial }{\partial t}}-{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)T=-\mathit{v}\cdot {\mathrm{\nabla }}_{\mathrm{h}}T+{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}^{\mathrm{2}}T;\begin{array}{l}\mathit{v}={\mathit{v}}_{\mathrm{h}}-{\mathit{v}}_{\mathrm{d}}\\ {\mathit{v}}_{\mathrm{d}}=\frac{\mathrm{1}}{\mathit{\rho }c}{\mathrm{\nabla }}_{\mathrm{h}}\left({\mathit{\kappa }}_{\mathrm{h}}\mathit{\rho }c\right),\end{array}\end{array}$$

where I have defined an effective diffusion velocity v_d and effective advection velocity v. Eq. (19) must be solved with the boundary conditions in Eq. (12).

The roles of the various terms are clearer if the equation is nondimensionalized. For this, note that if we include the boundary conditions, the anomaly temperature is entirely determined by the dimensional quantities κ, s, ρ, and c. From these, there is a unique dimensional combination τ(x) with dimensions of time; we will see that this controls the relaxation of the system back to energy balance, and it is a “relaxation time” (for the zero-dimensional model, energy balance is the same as thermodynamic equilibrium). Using κ_v yields

$$\begin{array}{}\text{(20)}& \mathit{\tau }={\mathit{\kappa }}_{\mathrm{v}}{\left(\mathit{\rho }cs\right)}^{\mathrm{2}};l_{\mathrm{v}}={\left(\mathit{\tau }{\mathit{\kappa }}_{\mathrm{v}}\right)}^{\mathrm{1}/\mathrm{2}}={\mathit{\kappa }}_{\mathrm{v}}\mathit{\rho }cs,\end{array}$$

where l_v(x) is the vertical relaxation length of the surface energy balance processes. In the next section, I give some rough parameter estimates. We may also define the horizontal diffusion length l_h, speed V, nondimensional (square root) diffusivity ratio β, and nondimensional advection vector α:

$$\begin{array}{}\text{(21)}& \mathit{\alpha }={\displaystyle \frac{\mathit{v}}{V}};V={\displaystyle \frac{l_{\mathrm{h}}}{\mathit{\tau }}};l_{\mathrm{h}}=(\mathit{\tau }{\mathit{\kappa }}_{\mathrm{h}}{)}^{\mathrm{1}/\mathrm{2}}=\mathit{\beta }{\mathit{\kappa }}_{\mathrm{h}}\mathit{\rho }cs;\mathit{\beta }={\left({\displaystyle \frac{{\mathit{\kappa }}_v}{{\mathit{\kappa }}_{\mathrm{h}}}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

The continuity equation for energy becomes $\mathrm{\nabla }\cdot \left(\frac{\mathit{\beta }}{s}\mathit{\alpha }\right)=\mathrm{0}$. For global (zero-dimensional) models, τ has been estimated as 2.4–7.0 years (90 % confidence; Procyk et al., 2020), which is comparable to the classical exponential relaxation timescales mentioned above (Hebert, 2017), and in Sect. 3.3 it is estimated as τ≈ 2.75 years.

In order to understand the classical origin of fractional derivatives, it is helpful to consider the homogeneous Sellers-type (diffusive transport) heat equation, where τ, l_v, and l_h are constants and can thus be used to nondimensionalize the operators. The nondimensional t is therefore in terms of relaxation times, the nondimensional x in terms of diffusion lengths l_h, and the nondimensional z in terms of diffusion depths l_v. Taking s= 1 effectively uses a forcing F with dimensions of temperature. The result is an equation with nondimensional operators acting on temperatures. In Part 1, I consider only the zero-dimensional equation, where “zero” refers to the number of horizontal dimensions (i.e., only vertical z and time t).

Using the dimensional parameters in Eqs. (20) and (21), we can write the equations as

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathit{Q}}_{\mathrm{h}}=-{\displaystyle \frac{l_{\mathrm{h}}}{s}}{\mathrm{\nabla }}_{\mathrm{h}}T+{\displaystyle \frac{\mathit{\alpha }}{s}}\left(T-T_{\mathrm{0}}\right);Q_z=-{\displaystyle \frac{l_{\mathrm{v}}}{s}}{\displaystyle \frac{\partial T}{\partial z}},\text{(23)}& {\displaystyle }\mathit{\tau }{\displaystyle \frac{\partial T}{\partial t}}=-\mathit{\zeta }T-l_{\mathrm{v}}s{\displaystyle \frac{\partial Q_z}{\partial z}};\begin{array}{l}\mathit{\zeta }=l_{\mathrm{h}}s{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(\frac{\mathit{\alpha }-l_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}}{s}\right)\\ \mathit{\zeta }T=l_{\mathrm{h}}s{\mathrm{\nabla }}_{\mathrm{h}}\cdot {\mathit{Q}}_{\mathrm{h}},\end{array}\end{array}$$

where ζ is the dimensionless horizontal transport operator. The reference temperature T₀ was ignored, justified by either taking it to be zero or assuming ${\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(s^{-\mathrm{1}}\mathit{\alpha }\right)=\mathrm{0}$, which is true if β is constant.

If the advection is chosen appropriately (as in Eq. 24 below), then we may write the horizontal transport operator in the form

$$\begin{array}{}\text{(24)}& \mathit{\zeta }=-s{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left({\displaystyle \frac{l_{\mathrm{h}}}{s}}\right){\mathrm{\nabla }}_{\mathrm{h}};\mathit{\alpha }=s{\mathrm{\nabla }}_{\mathrm{h}}\left({\displaystyle \frac{l_{\mathrm{h}}}{s}}.\right)\end{array}$$

This is convenient for comparing the HEBE with the 1D B-S equations on a sphere in Part 2 (Sect. 2.3), and it avoids the unphysical negative diffusivities reported by Sellers. Following North and Kim (2017), in spherical geometry, introduce D_F, which is the diffusion constant per radian:

$$\begin{array}{}\text{(25)}& \mathit{\zeta }=-sR{\mathrm{\nabla }}_{\mathrm{h}}\cdot D_F{\mathrm{\nabla }}_{\mathrm{h}};D_F={\displaystyle \frac{l_{\mathrm{h}}\left(\mathit{x}\right)}{Rs\left(\mathit{x}\right)}}={\mathit{\kappa }}_{\mathrm{h}}{\displaystyle \frac{\mathit{\beta }\mathit{\rho }c}{R}},\end{array}$$

where R is the Earth's radius.

3.1.2 Parameter estimates

Before proceeding, it is useful to get a feel for typical values of the parameters in the equations. In Sect. 2.3 and Appendix A of Part 2, I combine these parameter estimates with analyses of monthly space–time temperature anomalies in order to analyze which terms in the equations are dominant at different timescales; the following are order-of-magnitude estimates. The basic parameters are the horizontal diffusivity κ_h, the volumetric specific heat ρc, the sensitivity s, vertical diffusivity κ_v, κ_h, and relaxation time τ. They can be estimated as follows.

• a.

  Volumetric specific heat ρc. Ocean and land values are similar; the values for water and soil are $\mathit{\rho }c\approx \mathrm{4}\times {\mathrm{10}}^{\mathrm{6}}$ and ≈ 1 × 10⁶ J $/$ (m³ K), respectively. The soil value depends on moisture and soil type; this is an order-of-magnitude estimate.

• b.

  Climate sensitivity s. Using the CO₂ doubling value 3 ± 1.5 K, a 90 % confidence interval, and 3.71 W m² for CO₂ doubling, the global mean value is s≈ 0.8 ± 0.4 K W m², with regional values a factor of ≈ 2 higher or lower (IPCC AR5), yielding ρcs≈ 3 × 10⁶ s/m.

• c.

  Relaxation time τ. Based on responses to anthropogenic forcings since 1880, Hebert (2017), Hébert et al. (2021), and Procyk et al. (2020) give the global estimate τ ≈ 10⁸ s (≈ 4 years). This is comparable to the relaxation times for global box models.

• d.

  Horizontal diffusivity κ_h. As detailed in Part 2 (Sect. 2.3), North et al. (1981) and North and Kim (2017) use a diffusion constant per radian D_F (Eq. 25) combined with global-scale climatological forcing and temperature data to estimate a global thermal conductivity of K= 4.1 × 10⁶ Wm⁻¹/K from which we estimate the horizontal (eddy) diffusivity as κ_h=K /(ρc) ≈ 1 m²/s. Sellers (1969) gives values about 100 times larger for the ocean.

• e.

  Vertical diffusivity κ_v. The vertical diffusivity is not used in the usual energy balance models; however, in climate models, ocean values of κ_v≈ 10⁻⁴ m²/s are typical (Houghton et al., 2001). For soil, rough values are κ_v≈ 10⁻⁶ m²/s (wet) and κ_v ≈ 10⁻⁷ m²/s (dry); see Márquez et al. (2016). Alternatively, we can use ${\mathit{\kappa }}_{\mathrm{v}}=\mathit{\tau }/(\mathit{\rho }cs{)}^{\mathrm{2}}$ and the global estimates of τ≈ 10⁸ s to obtain κ_v ≈ 10⁻⁵ m²/s, which is close to the model values.

• f.

  Diffusion depth l_v. Using l_v=κ_vρcs, we find for the ocean and soils l_v≈ 300 m and ≈ 3–10 m, respectively. Using the global estimates, κ_v≈ 10⁻⁵–10⁻⁴ m²/s yields l_v≈ 30–100 m.

• g.

  Diffusion length l_h. This is a key parameter: $l_{\mathrm{h}}={\left(\mathit{\tau }{\mathit{\kappa }}_{\mathrm{h}}\right)}^{\mathrm{1}/\mathrm{2}}=\mathit{\beta }{\mathit{\kappa }}_{\mathrm{h}}\mathit{\rho }cs$ (Eq. 21). Using $l_{\mathrm{h}}={\left({\mathit{\kappa }}_{\mathrm{h}}{\mathit{\kappa }}_{\mathrm{v}}\right)}^{\mathrm{1}/\mathrm{2}}\mathit{\rho }cs$, l_h≈ 30 km (ocean) and 3 km (land). Using $l_{\mathrm{h}}={\left(\mathit{\tau }{\mathit{\kappa }}_{\mathrm{h}}\right)}^{\mathrm{1}/\mathrm{2}}$and κ_h≈1 m²/s yields a global estimate of l_h≈ 10 km.

• h.

  Diffusive-based velocity parameter V. $V\approx l_{\mathrm{h}}/\mathit{\tau }\approx$ 3 × 10⁻³–3 × 10⁻⁴ m/s.

• i.

  Nondimensional advection velocity α. The best transport model – diffusive, advective, or both – is not clear; therefore, let us estimate the magnitude of the advective velocity v assuming that it dominates the transport. The appropriate value is not obvious since most models just use eddy diffusivity – not advection – for transport. One way, for example (Warren and Schneider, 1979), is to note that typical meridional heat fluxes are of the order of 100 W m² over meridional bands whose temperature gradients ΔT are several Kelvin. If this heat is transported by advection, it implies $v\approx Q_{\mathrm{a}}/\left(\mathit{\rho }c\mathrm{\Delta }T\right)\approx$ 10⁻⁵–10⁻⁴ m/s (Eq. 4); hence, using V≈ 10⁻⁴ m/s (above), we find $\mathit{\alpha }=v/V\approx$ 0.1–1.

3.1.3 The nondimensional equations

With z and t in dimensionless form, the homogeneous zero-dimensional heat equation is

$$\begin{array}{}\text{(26)}& \left({\displaystyle \frac{\partial }{\partial t}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)T\left(t;z\right)=\mathrm{0}.\end{array}$$

I use the following notation: the first argument is t, then (when applicable) horizontal space, then a semicolon followed by the depth z. The transfer is confined to the semi-infinite region z≤ 0 with boundary conditions $T\left(t;-\mathrm{\infty }\right)=\mathrm{0}$ (bottom). The system is forced by the conductive–radiative surface boundary condition at z= 0 (the top):

$$\begin{array}{}\text{(27)}& {\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_{z=\mathrm{0}}+T\left(t;\mathrm{0}\right)=F\left(t\right).\end{array}$$

For initial conditions, in this section, the forcing is “turned on” at t > 0 (i.e., T(t; z)= 0 for t≤ 0), allowing use of Laplace transforms (see Sect. 3.3 for Fourier methods).

Performing a Laplace transform (LT) of the heat equation, we obtain

$$\begin{array}{}\text{(28)}& \left({\displaystyle \frac{d^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}}-p\right)\widehat{T}\left(p;z\right)=-T\left(\mathrm{0};z\right)=\mathrm{0},\end{array}$$

where the circumflex indicates the Laplace transform in time (with conjugate variable p):

$$\begin{array}{}\text{(29)}& \widehat{T}(p;z)=A\left(p\right)e^{\sqrt{p}z}+B\left(p\right)e^{-\sqrt{p}z},\end{array}$$

where A and B are determined by the BCs. Since the temperature at depth (z≪0) remains finite, we must have B= 0, and hence

$$\begin{array}{}\text{(30)}& \widehat{T}(p;z)=A\left(p\right)e^{\sqrt{p}z}.\end{array}$$

To determine A(p), Laplace-transform the surface boundary condition:

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{d}\widehat{T}}{\mathrm{d}z}}{|}_{z=\mathrm{0}}+\widehat{T}(p;\mathrm{0})=\widehat{F}\left(p\right);F\left(t\right)\stackrel{\mathrm{LT}}{\leftrightarrow }\widehat{F}\left(p\right),\end{array}$$

yielding

$$\begin{array}{}\text{(32)}& A\left(p\right)={\displaystyle \frac{\widehat{F}\left(p\right)}{\mathrm{1}+\sqrt{p}}}.\end{array}$$

It is more convenient to determine the response G_δ(t;z) to the impulse forcing F(t)=δ(t), which is the impulse Green's function. Using Eqs. (30) and (32) we obtain

$$\begin{array}{}\text{(33)}& \widehat{G_{\mathit{\delta }}}\left(p;z\right)={\displaystyle \frac{e^{\sqrt{p}z}}{\mathrm{1}+\sqrt{p}}};F\left(t\right)=\mathit{\delta }\left(t\right)\stackrel{\mathrm{LT}}{\leftrightarrow }\widehat{F}\left(p\right)=\mathrm{1}.\end{array}$$

The above assumes that the subsurface is infinitely deep. If instead it has a finite thickness L, and we take the bottom boundary condition as $T\left(t;-L\right)=\mathrm{0}$ (rather than $T\left(t;-\mathrm{\infty }\right)=\mathrm{0})$, then $B\left(p\right)\approx O\left(e^{-\mathrm{2}L\sqrt{p}}\right)$ and $\widehat{G_{\mathit{\delta }}}\left(p;\mathrm{0}\right)=\frac{\mathrm{1}}{\mathrm{1}+\sqrt{p}}-\frac{\mathrm{2}{e}^{-\mathrm{2}L\sqrt{p}}\sqrt{p}}{{\left(\mathrm{1}+\sqrt{p}\right)}^{\mathrm{2}}}+O\left(e^{-\mathrm{4}L\sqrt{p}}\right)$ so that the influence of the bottom condition on the surface decreases exponentially fast as its depth L increases. Physically, as long as the depth is of the order of a few diffusion depths (estimated as ≈ 100 m in the ocean, ≈ 10 m for land), the semi-infinite geometry assumption is unimportant. In the following, I therefore ignore any finite thickness corrections.

Taking the inverse Laplace transform of Eq. (33) we obtain the integral representation:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{G}_{\mathit{\delta }}\left(t;z\right)& ={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathit{\zeta }{e}^{-{\mathit{\zeta }}^{\mathrm{2}}t}}{\mathrm{1}+{\mathit{\zeta }}^{\mathrm{2}}}}\left(-\sin z\mathit{\zeta }+\mathit{\zeta }\cos z\mathit{\zeta }\right)\mathrm{d}\mathit{\zeta }\\ & \stackrel{\mathrm{LT}}{\leftrightarrow }\widehat{G_{\mathit{\delta }}}\left(p;z\right)={\displaystyle \frac{e^{\sqrt{p}z}}{\mathrm{1}+\sqrt{p}}},\end{array}\end{array}$$

with z≤0, where I have used contour integration on the Bromwich integral.

3.1.4 The surface temperature

For the surface, the integral (Eq. 34) can be expressed with the help of higher mathematical functions.

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{G}_{\mathrm{0},\mathrm{1}/\mathrm{2}}& \left(t;\mathrm{0}\right)=G_{\mathit{\delta }}\left(t;\mathrm{0}\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathit{\pi }t}}}-e^t\mathrm{erfc}\sqrt{t}\\ & \stackrel{\mathrm{LT}}{\leftrightarrow }{\widehat{G}}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(p;\mathrm{0}\right)=\widehat{G_{\mathit{\delta }}}\left(p;\mathrm{0}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\sqrt{p}}};\\ & \mathrm{erfc}\left(z\right)={\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}}\underset{z}{\overset{\mathrm{\infty }}{\int }}{e}^{-u^{\mathrm{2}}}\mathrm{d}u\end{array}\end{array}$$

G_δ(t;0) is the h= $\mathrm{1}/\mathrm{2}$ impulse response Green's function, also denoted as $G_{\mathrm{0}}{,}_{\mathrm{1}/\mathrm{2}}$; the 0 is for the 0th integral of the impulse, and $\mathrm{1}/\mathrm{2}$ is for the order of the derivative for its equation; see below. It is sometimes called a “generalized exponential”, which is itself expressed in terms of Mittag–Leffler functions.

For long times after an impulse, the response $G_{\mathit{\delta }}\left(t;\mathrm{0}\right)\approx t^{-\mathrm{3}/\mathrm{2}}$ (t≫1, Eq. 37 below) so that the system rapidly returns to its original temperature. It is more interesting to consider the response of the system to a step (Heaviside) forcing F(t)=Θ(t) (= 1, for t > 0, = 0 for t≤0) after which the system eventually attains a new energy balance (for the zero-dimensional model, this corresponds to thermodynamic equilibrium). Since $\mathrm{\Theta }\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\delta }\left(u\right)\mathrm{d}u$, we have the step response $G_{\mathrm{\Theta }}\left(t;z\right)=\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathit{\delta }}\left(u;z\right)\mathrm{d}u$ (also denoted $G_{\mathrm{1}}{,}_{\mathrm{1}/\mathrm{2}}$ in Eq. 36), and $G_{\mathrm{\Theta }}\left(t;\mathrm{0}\right)\approx \mathrm{1}-\frac{\mathrm{1}}{\sqrt{\mathit{\pi }t}}$ (Eq. 37), i.e., a slow power-law approach to thermodynamic equilibrium. Figs. 3 and 4 show this at different times and depths. With unit step forcing, the boundary condition (Eq. 27) indicates that the fraction of the heat flux that is transformed into longwave radiation is equal to the temperature with unit forcing. Therefore, the z= 0 curve in Fig. 3 shows that at first, all the forcing flux is conducted into the subsurface, but this fraction rapidly vanishes as the surface approaches equilibrium. At equilibrium, the temperature has increased so that the shortwave and longwave fluxes are once again in balance and there is no longer any conductive flux.

Figure 3 The nondimensional temperature as a function of nondimensional time for various nondimensional depths with a step forcing: G_Θ(t; z) (obtained by integrating Eq. 34 in time). The (top) surface curve can be interpreted as the fraction of the forcing that is conductive. At first, all the forcing is conductive with no radiation; eventually, all the fluxes are radiative, the system reaches a new thermodynamic equilibrium, and there is no conductive heat flux.

Figure 4 Contours of nondimensional temperature as a function of nondimensional time and depth after a step function forcing (G_Θ(t; z)).

For future reference, I give the corresponding step response $G_{\mathrm{1},\mathrm{1}/\mathrm{2}}$ =G_Θ, which is the integral of $G_{\mathrm{0},\mathrm{1}/\mathrm{2}}$ that describes relaxation to energy balance (for this model, thermodynamic equilibrium) when F is a step function. Similarly, the ramp (linear forcing) response $G_{\mathrm{2},\mathrm{1}/\mathrm{2}}$ is the integral of the step response, the second integral of the Dirac.

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{G}_{\mathrm{1},\mathrm{1}/\mathrm{2}}\left(t\right)& =G_{\mathrm{\Theta },\mathrm{1}/\mathrm{2}}\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(s\right)\mathrm{d}s\\ & =\mathrm{1}-e^t\mathrm{erfc}\left(t^{\mathrm{1}/\mathrm{2}}\right)\end{array}\end{array}$$

$$\begin{array}{rl}{G}_{\mathrm{2},\mathrm{1}/\mathrm{2}}\left(t\right)& =\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathrm{1},\mathrm{1}/\mathrm{2}}\left(s\right)\mathrm{d}s\\ & =\mathrm{1}-\mathrm{2}\sqrt{{\displaystyle \frac{t}{\mathit{\pi }}}}+t-e^t\mathrm{erfc}\left(t^{\mathrm{1}/\mathrm{2}}\right)\end{array}$$

For small and large t, we have the following.

$$\begin{array}{}\text{(37)}& G_{\mathrm{0},\mathrm{1}/\mathrm{2}}\left(t\right)=G_{\mathit{\delta },\mathrm{1}/\mathrm{2}}\left(t\right)\approx \begin{array}{l}\frac{\mathrm{1}}{\sqrt{\mathit{\pi }t}}-\mathrm{1}+\mathrm{2}\sqrt{\frac{t}{\mathit{\pi }}}-t\\ +\frac{\mathrm{4}}{\mathrm{3}}t\sqrt{\frac{t}{\mathit{\pi }}}-t-\mathrm{\dots }\ll \mathrm{1}\\ \frac{\mathrm{1}}{\mathrm{2}t\sqrt{\mathit{\pi }t}}-\frac{\mathrm{3}}{\mathrm{4}}\frac{\mathrm{1}}{t^{\mathrm{2}}\sqrt{\mathit{\pi }t}}+t+\mathrm{\dots }\gg \mathrm{1}\end{array}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{G}_{\mathrm{1},\mathrm{1}/\mathrm{2}}\left(t\right)=G_{\mathrm{\Theta },\mathrm{1}/\mathrm{2}}\left(t\right)\approx \begin{array}{l}\mathrm{2}\sqrt{\frac{t}{\mathit{\pi }}}-t+\frac{\mathrm{4}}{\mathrm{3}}\frac{t^{\mathrm{3}/\mathrm{2}}}{\sqrt{\mathit{\pi }}}-t-\mathrm{\dots }\ll \mathrm{1}\\ \mathrm{1}-\frac{\mathrm{1}}{\sqrt{\mathit{\pi }t}}+\frac{\mathrm{1}}{\mathrm{2}t\sqrt{\mathit{\pi }t}}-t-\mathrm{\dots }\gg \mathrm{1}\end{array}\\ {\displaystyle }& {\displaystyle }{G}_{\mathrm{2},\mathrm{1}/\mathrm{2}}\left(t\right)\approx \begin{array}{l}\frac{\mathrm{4}}{\mathrm{3}}t\sqrt{\frac{t}{\mathit{\pi }}}-\frac{t^{\mathrm{2}}}{\mathrm{2}}+\frac{\mathrm{8}}{\mathrm{15}}{t}^{\mathrm{2}}\sqrt{\frac{t}{\mathit{\pi }}}-\frac{t^{\mathrm{3}}}{\mathrm{6}}+t+\mathrm{\dots }\ll \mathrm{1}\\ t+\mathrm{1}-\mathrm{2}\sqrt{\frac{t}{\mathit{\pi }}}-\frac{\mathrm{1}}{\sqrt{\mathit{\pi }t}}+\frac{\mathrm{1}}{\mathrm{2}t\sqrt{\mathit{\pi }t}}-t-\mathrm{\dots }\gg \mathrm{1}\end{array}\end{array}$$

The asymptotic equation for the step response ($G_{\mathrm{1},\mathrm{1}/\mathrm{2}}$) shows that equilibrium is approached slowly as $t^{-\mathrm{1}/\mathrm{2}}$. It is this power-law step response (empirically with $\approx t^{-\mathrm{0.5}})$ that was discovered semi-empirically by Hebert (2017) and Lovejoy et al. (2017, 2021) and was successfully used for climate projections through 2100 (Hébert et al., 2021). Similarly, $\approx t^{-\mathrm{0.4}}$ behavior was used for macroweather (monthly, seasonal) forecasts close to the short-time $t^{-\mathrm{1}/\mathrm{2}}$ impulse response expansion (Lovejoy et al., 2015; Del Rio Amador and Lovejoy, 2019).

If we take this as a model of the global temperature, we can use the ramp Green's function to estimate the ratio of the equilibrium climate response (ECS) to the transient climate response (TCR); we find TCR $/$ ECS$=G_{\mathrm{2},\mathrm{1}/\mathrm{2}}\left(\mathrm{\Delta }t\right)/\mathrm{\Delta }t$, where Δt is the nondimensional time over which (for the TCR) the linear forcing acts. Using τ= 4 years and the standard Δt= 70 years for the TCR ramp, we find the plausible ratio TCR/ECS ≈ 0.78.

3.1.5 Comparison with temperature forcing boundary conditions

It is interesting to compare this with the classical surface boundary condition when the system is forced by the surface temperature; an alternative – periodic surface heat forcing – is discussed in Sect. 3.3. If the surface (z= 0) boundary condition T_force(t) is imposed,

$$\begin{array}{}\text{(38)}& T_{\mathrm{temp}}\left(t;\mathrm{0}\right)=T_{\mathrm{force}}\left(t\right),\end{array}$$

then there will be vertical surface gradients that imply that heat is conducted through the surface. To obtain the impulse response Green's function, take T_force(t)=δ(t), and repeating the Laplace transform approach, we obtain A(p)= 1 (Eq. 31 with no derivative term). This yields the following Laplace transform pairs for the impulse and step Green's function.

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& G_{\mathrm{temp},\mathit{\delta }}\left(t;z\right)={\displaystyle \frac{z{e}^{-z^{\mathrm{2}}t}}{\mathrm{2}\sqrt{\mathit{\pi }{t}^{\mathrm{3}}}}}\stackrel{\mathrm{LT}}{\leftrightarrow }{\widehat{G}}_{\mathrm{temp},\mathit{\delta }}(p;z)=e^{\sqrt{p}z}\\ & G_{\mathrm{temp},\mathrm{\Theta }}(t;z)=\mathrm{1}+\mathrm{erf}\left({\displaystyle \frac{z}{\mathrm{2}\sqrt{t}}}\right)\\ & \stackrel{\mathrm{LT}}{\leftrightarrow }{\widehat{G}}_{\mathrm{temp},\mathrm{\Theta }}(p;z)={\displaystyle \frac{e^{\sqrt{p}z}}{p}}\end{array}\end{array}$$

In the context of the Earth's temperature, using heat conduction (not temperature) boundary conditions, Brunt (1932) obtained the analogous classical formula noting that “this solution is given in any textbook”.

These classical Green's functions provide useful comparisons with the conductive–radiative BCs. For example, integrating Eq. (34) with respect to time and simplifying, we obtain the following.

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& \mathrm{\Delta }{G}_{\mathrm{\Theta }}\left(t;z\right)=G_{\mathrm{\Theta },\mathrm{temp}}\left(t;z\right)-G_{\mathrm{\Theta }}\left(t;z\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{e^{-{\mathit{\zeta }}^{\mathrm{2}}t}{e}^{iz\mathit{\zeta }}\mathrm{d}\mathit{\zeta }}{\left(\mathrm{1}+i\mathit{\zeta }\right)}};\begin{array}{l}t\ge \mathrm{0}\\ z\le \mathrm{0}\end{array}\end{array}\end{array}$$

Since the step response G_Θ describes the approach to thermodynamic equilibrium, ΔG_Θ(t;z) (Fig. 5) succinctly expresses the differences between the temperature and conductive–radiative forced boundary conditions. The leading large t approximation to the integral in Eq. (40) is $\mathrm{\Delta }{G}_{\mathrm{\Theta }}\left(t;z\right)\approx e^{-\frac{z^{\mathrm{2}}}{\mathrm{4}t}}/\sqrt{\mathit{\pi }t}$ so that, although they both slowly approach each other and eventually attain equilibrium, the differences are important (especially in the diffusion layer, z≈ < 1) and they decay very slowly with time and depth; I discuss this further in Sect. 3.3.

Figure 5 The difference ΔG_Θ(t;z) between the classical (temperature-forced) and radiatively forced step response functions over the diffusion depth (nondimensional z= 0 to −1). The top shows the surface (z= 0), and the curves from the top to bottom are at depths z= 0, −0.1, −0.2, −0.3, … −1. While the difference is large over the relaxation time (up to nondimensional t= 1), we see that they both slowly converge to thermodynamic equilibrium at large t.

3.1.6 Surface temperatures, fractional derivatives, and the HEBE

Let us now introduce the hth-order fractional derivative ${}_{t_{\mathrm{0}}}{D}_t^h$ to represent the fractional derivative order h of an arbitrary function f over the domain from t₀ to t:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& {}_{t_{\mathrm{0}}}{D}_t^{h}f={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-h\right)}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\left(t-u\right)}^{-h}{f}^{\prime }\left(u\right)\mathrm{d}u;\\ & f^{\prime }\left(u\right)={\displaystyle \frac{\mathrm{d}f}{\mathrm{d}u}};\mathrm{0}\le h\le \mathrm{1}.\end{array}\end{array}$$

Fractional derivatives of order h are most commonly interpreted in the Caputo (as above) or Riemann–Liouville sense (Podlubny, 1999). For h≤1, the main case of interest here, the distinction is not important and they most commonly use t₀= 0. Fractional derivatives and their inverses, fractional integrals (with h < 0), are thus power-law-weighted convolutions; fractional integrals of noises are often associated with long-memory stochastic processes. Many studies have found long memories in macroweather (Blender and Fraedrich, 2003; Bunde et al., 2005; Rybski et al., 2006; Varotsos et al., 2013), and Gaussian-noise-forced models (fractional Gaussian noise) have been proposed as models of internally forced (macroweather) temperature variability (Rypdal and Rypdal, 2014; Lovejoy, 2015; Del Rio Amador and Lovejoy, 2019; Del Rio Amador and Lovejoy, 2021a).

Most applications of fractional derivatives are for forcings that start at $t=t_{\mathrm{0}}=$ 0 (i.e., F= 0 for t≤0) (see Miller and Ross, 1993, and Podlubny, 1999) and are convenient for deterministic forcings; however, they singularize t= 0, whereas we often wish to include periodic or statistically stationary internal stochastic forcings so that $F\left(-\mathrm{\infty }\right)=\mathrm{0}$ (or in the periodic case, the mean over a cycle = 0) is more convenient, in which case take $t_{\mathrm{0}}=-\mathrm{\infty }$ and hence $T_{\mathrm{s}}\left(t=-\mathrm{\infty }\right)=$ 0 (or periodic). As discussed in Lovejoy (2019b), this corresponds to the semi-infinite-range Weyl fractional derivative. Deterministic, stochastic, and periodic forcings can be combined into a single framework simply by using the Weyl derivatives with, for example, the deterministic part of the forcing starting at t= 0 (with the deterministic F(t)= 0 for t≤0) and the stochastic forcing at $t=-\mathrm{\infty }$. These fractional derivatives have the following transformation properties:

$$\begin{array}{}\text{(42)}& \begin{array}{rl}& {}_{\mathrm{0}}{D}_t^h\stackrel{\mathrm{LT}}{\leftrightarrow }{p}^h\\ & {}_{-\mathrm{\infty }}{D}_t^h\stackrel{\mathrm{FT}}{\leftrightarrow }{\left(i\mathit{\omega }\right)}^h,\end{array}\end{array}$$

where ω is the Fourier conjugate to t (see, e.g., Miller and Ross, 1993; Podlubny, 1999). In this Part 1 (except for Sect. 3.3), I consider deterministic forcings; putting t₀= 0 in Eq. (41) and using ${}_{\mathrm{0}}{D}_t^{\mathrm{1}/\mathrm{2}}\stackrel{\mathrm{LT}}{\leftrightarrow }\sqrt{p}$ (h= 1/2 in Eq. 42), we obtain the HEBE for the surface temperature Green's function.

$$\begin{array}{}\text{(43)}& \left({}_{\mathrm{0}}{D}_t^{\mathrm{1}/\mathrm{2}}+\mathrm{1}\right)G_{\mathit{\delta }}\left(t;\mathrm{0}\right)=\mathit{\delta }\left(t\right)\stackrel{\mathrm{LT}}{\leftrightarrow }\left(\sqrt{p}+\mathrm{1}\right){\widehat{G}}_{\mathit{\delta }}(p;\mathrm{0})=\mathrm{1}\end{array}$$

This proves that the surface temperatures implied by the heat equation with conductive–radiative boundary conditions can be determined directly from the HEBE using the same Green's function. For the dimensional equations, the surface temperature therefore satisfies the dimensional HEBE:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}& {\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}{}_{\mathrm{0}}{D}_t^{\mathrm{1}/\mathrm{2}}{T}_{\mathrm{s}}+T_{\mathrm{s}}=sF\left(t\right);\\ & T_{\mathrm{s}}\left(t\right)=s\underset{\mathrm{0}}{\overset{t}{\int }}{G}_{\mathit{\delta }}\left({\displaystyle \frac{t-u}{\mathit{\tau }}};\mathrm{0}\right)F\left(u\right){\displaystyle \frac{\mathrm{d}u}{\mathit{\tau }}},\end{array}\end{array}$$

where the surface temperature is $T_{\mathrm{s}}\left(t\right)=T(t;\mathrm{0})$.

This HEBE for the surface temperature could be regarded as a significant nonclassical example of the Mori–Zwanzig formalism (Gottwald et al., 2017; Mori, 1965; Zwanzig, 1973, 2001) and empirical model reduction formalisms (Ghil and Lucarini, 2020), whereby memory effects arise if we only look at one part of the system, ignoring the others. In the HEBE, the surface temperature is analogously expressed directly in terms of the forcing, ignoring the subsurface degrees of freedom. Although such memories are usually considered exponential and hence small, the HEBE shows that the classical continuum heat equation has, on the contrary, strong power-law memories. This points to serious limitations on conventional dynamical systems approaches to climate science that assume that the dynamical equations are integer-ordered with exponential memories. The HEBE shows that the (fundamental) radiatively exchanging components of the climate system will generally be characterized by long memories associated with fractional rather than integer-ordered derivatives. I develop this insight elsewhere.

3.2 The HEBE, zero-dimensional and box models, and Newton's law of cooling

Phenomenological models of the temperature based on the energy balance across a homogeneous surface may represent either the whole Earth or only a subregion. The former are global zero-dimensional energy balance models (sometimes called global energy balance models or GEBMs; see the review in McGuffie and Henderson-Sellers, 2005), whereas the latter may represent the balance across the surface of a homogeneous subsection: a “box”. The boxes have spatially uniform temperatures that store energy according to their heat capacity, density, and size. Often several boxes are used, mutually exchanging energy, and the basic idea can be extended to column models. Since the average Earth temperature can be modeled either as a single horizontally homogeneous box or by two or more vertically superposed boxes, in the following, “box model” refers to both global and regional models.

A key aspect of these models is the rate at which energy is stored and at which it is exchanged between the boxes. Stored heat energy is transferred across a surface, and it is generally postulated that its flux obeys Newton's law of cooling (NLC). The NLC is usually only a phenomenological model; it states that a body's rate of heat loss is directly proportional to the difference between its temperature and its environment. In these horizontally homogeneous models, it is only the heat energy per area (S) that is important so that the NLC can be written as

$$\begin{array}{}\text{(45)}& Q_{\mathrm{s}}={\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{Z}}\left(T_{\mathrm{eq}}-T\right).\end{array}$$

S is the heat in the body and Q_s is the heat flux across the surface into the body (see Fig. 6). T_eq is the equilibrium temperature, and Z is a transfer coefficient, which is the “thermal impedance” (units: m² K/W); its reciprocal Y is the surface “thermal admittance” (see the next section). Identifying the equilibrium temperature with T_eq(t)= sF(t) and using the dimensional surface boundary condition (Eq. 12), it is easy to check that a direct consequence of the HEBE's conductive–radiative boundary condition is that it also satisfies the NLC:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{Q}_{\mathrm{s},\mathrm{HEBE}}& ={\displaystyle \frac{\mathrm{d}{S}_{\mathrm{HEBE}}}{\mathrm{d}t}}=\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_{z=\mathrm{0}}={\displaystyle \frac{\left(T_{\mathrm{eq}}-T\right)}{s}};\\ & T_{\mathrm{eq}}=sF.\end{array}\end{array}$$

Unlike the usual phenomenological box applications that simply postulate the NLC, the HEBE satisfies it as a consequence of its energy-conserving surface boundary condition. Comparing Eqs. (41) and (42), we may also conclude that thermal impedance Z=s.

Figure 6 A schematic showing Newton's law of cooling (NLC) that relates the temperature difference across a surface to the heat flux crossing the surface, Q_s (into the surface). T_eq is the fixed outside temperature; heat will flow as long as the surface temperature T_s≠T_eq. Z is the thermal impedance (equal here to the climate sensitivity s). To apply the NLC, we need to relate the heat flux to the surface temperature. The lower left shows the consequence of applying the heat equation with conductive–radiative BCs, and the lower right shows the phenomenological assumption made by box models. The arrows represent heat fluxes, hence the factor s in the denominators. The system is assumed to be horizontally homogeneous with a subsurface much thicker than the diffusion depth.

While the HEBE and box models both obey the NLC, the relationships between the surface heat flux Q_s= dS $/$ dt and the surface temperature T are quite different. For example, for forcings starting at time t=t₀, using the HEBE we have

$$\begin{array}{}\text{(47)}& \begin{array}{r}{Q}_{\mathrm{s},\mathrm{HEBE}}={\displaystyle \frac{\mathrm{d}{S}_{\mathrm{HEBE}}}{\mathrm{d}t}}={\displaystyle \frac{{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}}{s}}{}_{t_{\mathrm{0}}}{D}_t^{\mathrm{1}/\mathrm{2}}T;\\ \mathit{\tau }=\mathit{\rho }cs{l}_{\mathrm{v}};l_{\mathrm{v}}={\mathit{\kappa }}_{\mathrm{v}}\mathit{\rho }cs.\end{array}\end{array}$$

Although this relation between surface heat fluxes and temperatures has been known for some time (Babenko, 1986; Podlubny, 1999; see, e.g., Sierociuk et al., 2013, 2015, for applications), to my knowledge, it has never been applied to conduction–radiative models, nor has it been combined with the NLC to yield the homogeneous HEBE. In comparison, box models satisfy

$$\begin{array}{}\text{(48)}& Q_{\mathrm{s},\mathrm{box}}={\displaystyle \frac{\mathrm{d}{S}_{\mathrm{box}}}{\mathrm{d}t}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{box}}}{s}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}};{\mathit{\tau }}_{\mathrm{box}}=\mathit{\rho }csL;L={\displaystyle \frac{C}{\mathit{\rho }c}},\end{array}$$

where L is the effective thickness of the surface layer, C is the specific heat per area, and τ_box is the classical EBE relaxation time. Geoffroy et al. (2013) used a two-box model to fit outputs of a dozen general circulation models (GCMs) and found τ_box≈ 4.1 ± 1.1 years (the mean and spread of 12 models) and ≈ 40–800 years for the second box, whereas the IPCC (2013) recommends a two-box model with relaxation scales of τ_box= 8.4 and 409 years. With the FEBE, Procyk et al. (2020) find h= 0.38 ± 0.05 and τ= 4.7 ± 2.3 years.

The HEBE and box heat transfer models can conveniently be compared and contrasted by placing them both in a more general common framework. Define the hth-order heat storage as

$$\begin{array}{}\text{(49)}& S_{\mathrm{h}}\left(t\right)={\displaystyle \frac{{\mathit{\tau }}^h}{s\mathrm{\Gamma }\left(\mathrm{1}-h\right)}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}T\left(u\right){\left(t-u\right)}^{-h}\mathrm{d}u;\mathrm{0}\le h\le \mathrm{1}.\end{array}$$

If T(t₀)= 0 (this is equivalent to fixing the reference of the anomalies), then integrating by parts yields

$$\begin{array}{}\text{(50)}& S_{\mathrm{h}}\left(t\right)={\displaystyle \frac{{\mathit{\tau }}^{h-\mathrm{1}}}{s\mathrm{\Gamma }\left(\mathrm{1}-h\right)}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{T}^{\prime }\left(u\right){\left(t-u\right)}^{\mathrm{1}-h}\mathrm{d}u;\mathrm{0}\le h\le \mathrm{1}.\end{array}$$

Putting h= 1 yields the simple $S_{\mathrm{1}}\left(t\right)=T\left(t\right)/s$ so that S₁=S_box.

Over the interval t₀ to t, the fractional derivative of order h is defined as the ordinary derivative of the 1−h-order (Riemann–Liouville) fractional integral:

$$\begin{array}{}\text{(51)}& \begin{array}{rl}{}_{t_{\mathrm{0}}}{D}_t^{h}T& ={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({}_{t_{\mathrm{0}}}{D}_t^{h-\mathrm{1}}T\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}-h\right)}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\left(t-u\right)}^{-h}T\left(u\right)\mathrm{d}u\right];\\ & \mathrm{0}\le h\le \mathrm{1}.\end{array}\end{array}$$

Therefore, $S_{\mathrm{1}/\mathrm{2}}=S_{\mathrm{box}}$ and

$$\begin{array}{}\text{(52)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{h}}}{\mathrm{d}t}}=s^{-\mathrm{1}}{\mathit{\tau }}^h{}_{t_{\mathrm{0}}}{D}_t^{h}T;\begin{array}{l}{h}_{\mathrm{HEBE}}=\mathrm{1}/\mathrm{2};{\mathit{\tau }}_{\mathrm{FEBE}}=l_{\mathrm{v}}\mathit{\rho }cs,\\ h_{\mathrm{box}}=\mathrm{1};{\mathit{\tau }}_{\mathrm{box}}=L\mathit{\rho }cs.\end{array}\end{array}$$

Combining this with the NLC, in both cases we obtain

$$\begin{array}{}\text{(53)}& {\mathit{\tau }}^h{}_{t_{\mathrm{0}}}{D}_t^{h}T+T=sF.\end{array}$$

Hence, the box and HEBE models are special cases of the fractional-order energy balance equation (FEBE; Lovejoy, 2019a, b) derived phenomenologically in Lovejoy et al. (2021). Whereas the box model changes its heat content instantaneously with its current temperature (T(t)), at any moment, the energy stored in the HEBE model depends on the past temperatures, and since their weights fall off slowly – there is a long memory – it potentially depends on the temperature and hence energy stored in the distant past. Box or column models all have surfaces that exchanges heat both radiatively and conductively so that – contrary to standard practice – these surfaces should instead exchange heat fractionally with h= $\mathrm{1}/\mathrm{2}$ not h= 1. Note that when box interfaces with purely conductive heat exchanges are considered (without radiative transfer, e.g., between a “deep ocean” and “mixed layer” in global two-box model), then the thermal contact conductance that characterizes the interface is needed.

At a theoretical level, the advantage of the HEBE is that unlike, the box models, it is a direct consequence of the standard (energy-conserving) continuum heat equation combined with standard energy-conserving surface boundary conditions. It is therefore natural to ask if the h= 1 heat transfer (i.e., d$S_{\mathrm{1}}/\mathrm{d}t=$  ($C/s$)dT $/$ dt) can be derived from the heat transport equation. Returning to the nondimensional boundary condition $\left({\left.\frac{\partial T}{\partial z}\right|}_{z=\mathrm{0}}+T\left(t;\mathrm{0}\right)=F\left(t\right)\right)$, it is easy to verify that in order to recover h= 1 heat transfer, one must instead use ${\left.\frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}\right|}_{z=\mathrm{0}}+T\left(t;\mathrm{0}\right)=F\left(t\right)$. I therefore conclude that box model h= 1 transfer is not simultaneously compatible with the heat equation and energy balance boundary conditions.

To summarize, we are currently in the unsatisfactory position of having zero- and one-dimensional (box and Budyko–Sellers) energy balance equations, neither of which satisfy the correct radiative–conductive surface boundary conditions. For the box models, the consequence is that the energy storage processes have rapid (exponential) rather than slow (power-law) relaxation. For the Budyko–Sellers models, the consequence is that, at best, they are 1D, and even with this restriction, their time-dependent versions have derivatives of the wrong order (Part 2, Sect. 2.3). In comparison, the zero-dimensional HEBE is a consequence of correcting the Budyko–Sellers boundary conditions. It satisfies the NLC and corrects the order h by reducing it from the phenomenological value of h= 1 to h= $\mathrm{1}/\mathrm{2}$. As a bonus, in Part 2 we see that the HEBE can easily be extended from zero to two spatial dimensions, enlarging the scope of energy balance models while simultaneously eliminating these weaknesses.

3.3 Thermal impedance, complex climate sensitivities, and the annual cycle

3.3.1 Conductive versus conductive–radiative boundary conditions

Up until now, I have discussed forcing that is turned on at t= 0; this allowed for convenient solutions using Laplace transform methods. However, for forcing that is periodic or that is a stationary noise (i.e., the internal variability), Fourier techniques are more useful.

The first applications of Fourier techniques to the problem of radiative and conductive heat transfer into the Earth was by Brunt (1932) and Jaeger and Johnson (1953), who considered the (weather regime) diurnal cycle. I already mentioned that Brunt (1932) also considered step function heat forcing, which he claimed might be a plausible model of the diurnal cycle near sunset or sunrise. However, in zero-dimensional models, the long-time temperatures after step heat flux forcings are divergent (but not in 2D models; see Part 2) so that later in his paper Brunt considered periodic diurnal heat flux forcing with no net heat flux across the surface and used Fourier methods instead. In this classical diurnally forced problem, the periodic temperature response lags the forcing by a phase shift of $\mathit{\pi }/\mathrm{4}$ = 3 h. If we apply the same shift to the annual cycle – assuming that the Earth is forced by heat flux into its subsurface – the corresponding lag is 1.5 months ≈ 46 d, which is generally (at least for land) too long (we shall see that it corresponds to an infinite relaxation time).

Following Brunt (1932) and Jaeger and Johnson (1953), consider the response to a single Fourier component forcing (this is equivalent to Fourier analysis of the equation). In this case, assuming a periodic temperature response and substituting this into the 1D heat equation (time and depth, i.e., the dimensional version of Eq. 22), we find that the variation of amplitude with depth is

$$\begin{array}{}\text{(54)}& T\left(t;z\right)=T_{\mathrm{s}}{e}^{i\mathit{\omega }t}{e}^{\sqrt{\frac{i\mathit{\omega }}{{\mathit{\kappa }}_{\mathrm{v}}}}z};z\le \mathrm{0},\end{array}$$

where T_s is the amplitude of the surface temperature oscillations; it depends on the nature of the forcing, here on the boundary conditions (“s” is for surface). Following Brunt, using the classical heat surface heat forcing F_se^iωt as the surface boundary condition (with this forcing, F_s=Q_s is the heat crossing the surface and entering the system in the downward direction; see Figs. 1 and 6), we find

$$\begin{array}{}\text{(55)}& {\left.\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{heat}}}{\partial z}}\right|}_{z=\mathrm{0}}=F_{\mathrm{s}}{e}^{i\mathit{\omega }t},\end{array}$$

where “heat” is for heat forcing, and we obtain

$$\begin{array}{}\text{(56)}& T_{\mathrm{s},\mathrm{heat}}={\displaystyle \frac{F_{\mathrm{s}}}{\sqrt{i\mathit{\omega }{\left(\mathit{\rho }c\right)}^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{v}}}}}=Z\left(\mathit{\omega }\right)F_{\mathrm{s}};Z\left(\mathit{\omega }\right)={\displaystyle \frac{s}{\sqrt{i\mathit{\omega }\mathit{\tau }}}},\end{array}$$

where Z(ω) is the complex frequency-dependent thermal impedance, the reciprocal of the thermal admittance. For a given surface heat flux, Z(ω) quantifies the surface temperature response (I have written the impedance with the help of s in order to nondimensionalize the denominator). Thermal impedance and admittance are standard in areas of heat transfer engineering and were introduced into the problem of diurnal Earth heating by Byrne and Davis (1980). From Z(ω), we can thus easily understand the key (Brunt, 1932; Jaeger and Johnson, 1953) result that arg(Z(ω)) = arg($i^{-\mathrm{1}/\mathrm{2}})=-\mathit{\pi }$/4 (“arg” indicates the phase).

So far, this approach has only been applied to weather scales (the diurnal cycle). Let us now apply the same approach but with an eye to longer macroweather timescales, notably the annual cycle. The climate sensitivity is an emergent macroweather quantity determined by numerous feedbacks that over weather scales are quite nonlinear but over macroweather scales are considerably averaged (and, at least for GCMs, are already fairly linear; Hébert and Lovejoy, 2018). In any event, for the annual cycle we use radiative–conductive boundary conditions rather than the pure conductive ones used by Brunt.

Using conductive–radiative surface BCs with external forcing F_se^iωt yields

$$\begin{array}{}\text{(57)}& \begin{array}{rl}& F\left(t\right)=F_{\mathrm{s}}{e}^{i\mathit{\omega }t},\\ & F_{\mathrm{s}}=Q_{\mathrm{s}}+Q_{\mathrm{s},\mathrm{rad}}=s^{-\mathrm{1}}\left(\mathrm{1}+{\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}\right)T_{\mathrm{s}},\\ & Q_{\mathrm{s},\mathrm{rad}}=s^{-\mathrm{1}}{T}_{\mathrm{s}},\\ & {\left.Q_{\mathrm{s}}=\mathit{\rho }c{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial T}{\partial z}}\right|}_{z=\mathrm{0}}={\mathrm{s}}^{-\mathrm{1}}{\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}{T}_{\mathrm{s}},\end{array}\end{array}$$

where F_s is the radiative (downward) forcing radiative flux, and Q_s and Q_s,rad are the surface conductive (into the subsurface) and longwave radiative emission (away from the surface) fluxes, respectively. Solving, we obtain the same depth dependence (Eq. 54), but with the amplitude of the surface oscillations given by

$$\begin{array}{}\text{(58)}& T_{\mathrm{s}}=s\left(\mathit{\omega }\right)F_{\mathrm{s}};s\left(\mathit{\omega }\right)=Z\left(\mathit{\omega }\right)={\displaystyle \frac{s}{\mathrm{1}+{\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}}},\end{array}$$

where I have introduced the complex climate sensitivity s(ω), which by definition is equal to the complex thermal impedance Z(ω). In the context of the Earth's energy balance, it is more useful to think in terms of sensitivities than impedances so that below I use s(ω). With this, we obtain

$$\begin{array}{}\text{(59)}& Q_{\mathrm{s}}={\displaystyle \frac{s\left(\mathit{\omega }\right)}{s}}{\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}{F}_{\mathrm{s}};Q_{\mathrm{s},\mathrm{rad}}={\displaystyle \frac{s\left(\mathit{\omega }\right)}{s}}{F}_{\mathrm{s}}.\end{array}$$

Since arg($i^{\mathrm{1}/\mathrm{2}})=\mathit{\pi }$/4 (= 45^∘), we see that as mentioned earlier, the conductive and longwave radiative fluxes are out of phase by 45^∘, but the phase of the temperature lags the forcing by arg(s(ω)), which only reaches 45^∘ in the large τ limit (see Fig. 7).

Note that I could have deduced Eq. (59) directly by Fourier analysis of the HEBE using $\mathrm{FT}\left({}_{-\mathrm{\infty }}{D}_t^{\mathrm{1}/\mathrm{2}}\right)={\left(i\mathit{\omega }\right)}^{\mathrm{1}/\mathrm{2}}$, but the above allowed comparison with the results of the classical model. The Fourier method allows us to extend the complex climate sensitivity to the more general FEBE:

$$\begin{array}{}\text{(60)}& s_{\mathrm{h}}\left(\mathit{\omega }\right)={\displaystyle \frac{s}{\mathrm{1}+{\left(i\mathit{\omega }\mathit{\tau }\right)}^h}}.\end{array}$$

The usual EBE is the h= 1 special case.

Figure 7 The temperature phase lag (in months, the negative of argument of the complex climate sensitivity) using the complex climate sensitivity and annual cycle forcing (i.e., with ω= 2π rad per year) with τ in years. The line with short dashes (top) is the usual EBE (h= 1), the solid line is the (h= $\mathrm{1}/\mathrm{2}$) HEBE, and the line with long dashes is the classical heat forcing model, which is the large τ HEBE limit. All curves ignore any net horizontal heat transport. The data analyzed here yield τ ≈ 2.75 ± 0.8 years, but the actual phase is somewhat shorter due to horizontal heat transport.

Figure 8 Same as Fig. 7 except for the amplitude of the complex climate sensitivity to annual cycle forcing (i.e., with ω= 2π rad per year) with τ in years. The short dashed line (bottom) is the usual EBE (h = 1), the top line with long dashes is the classical heat forcing model, and the solid line is the (h= $\mathrm{1}/\mathrm{2}$) HEBE.

3.3.2 Empirical estimates of complex climate sensitivities

Figures 7 and 8 compare the phases and amplitudes of s(ω) for the HEBE with classical and conductive–convective boundary conditions (h= $\mathrm{1}/\mathrm{2}$) as well as the h= 1 EBE. The plots use ω= 2π rad per year. From Fig. 7, we see that, taking the empirical value τ≈ 4.7 years (Procyk et al., 2020), the HEBE lag is a little over a month. From the detailed maps in Donohoe et al. (2020) (see also Ziegler and Rehfeld, 2020) I estimate that in the extratropical regions over land, the summer temperature maximum is typically 30–40 d after the solstice but only 20–30 d after the maximum forcing (insolation). For ocean, it is 60–70 d after the solstice but only 30–40 d after the maximum insolation. The HEBE result is thus close to the observed lag between the summer solstice and maximum temperatures over most land areas.

In contrast, if the Brunt (1932) classical heat forcing result is used, we obtain π/4 = 1.5 months = 46 d, which is already too long for most of the globe, and the h= 1 EBE result (close to 3 months = 91 d) is much too long. Over the ocean, the lag is typically longer than over land, probably because of the strong albedo periodicity associated with seasonal ocean cloud cover (Stubenrauch et al., 2006; Donohoe et al., 2020). This delays the summer solstice forcing maximum over the ocean, potentially explaining the extra lag.

Although a complete analysis with modern data is out of our present scope, we can get a feel for the realism of this approach by using the zonally averaged (North and Coakley, 1979) Sellers model discussed in the review (North et al., 1981, updated in North et al., 1983) wherein most of the Earth follows the EBE phase lags of ≈ 90 d. The model uses a second-order Legendre polynomial to take into account the latitudinal variations and a sinusoidal annual cycle with empirically fit parameters that effectively zonally average over land and ocean. Empirical parameters are given for the albedo, top-of-the-atmosphere insolation, temperature, and outgoing IR radiation such that the global temperature maximum lags the solstice by 32.5 d (North and Coakley, 1979; North et al., 1983). An updated 2D version of the Sellers model has been used to estimate phase lags with respect to the solstice, finding lags of ≈ 90 d over oceans and ≈ 30–40 d over land (Zhuang et al., 2017; Ziegler and Rehfeld, 2020).

Before continuing, recall that the zero-dimensional theory discussed here assumes that radiative flux imbalances are all stored; it ignores the divergence of the horizontal heat transport, which according to Trenberth et al. (2009) is small even though the heat fluxes may be significant. Although, at least for temperature anomalies, I argue that this effect is mostly important at small scales, the magnitude of horizontal heat divergence at macroweather scales is not well known and is presumably quite variable from place to place depending on (inhomogeneous) local transport parameters (see Part 2). A simple way to parameterize the transport is to maintain the assumption that the Earth has homogeneous parameters and to assume that the transport is due to horizontally inhomogeneous forcing. In Part 2, I show that for a horizontal wavenumber k, the effect of horizontal transport is to modify the storage term as ${\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}\to {\left(i\mathit{\omega }\mathit{\tau }+{\left(l_{\mathrm{h}}k\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}$; therefore, for pure periodic horizontal forcing,

$$\begin{array}{}\text{(61)}& \begin{array}{rl}& Q_{\mathrm{s},\mathrm{h}}={\displaystyle \frac{s_{\mathrm{h}}\left(\mathit{\omega }\right){\left(i\mathit{\omega }\mathit{\tau }+{\left(l_{\mathrm{h}}k\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}}{s}}{F}_{\mathrm{s}};\\ & Q_{\mathrm{s},\mathrm{rad}}={\displaystyle \frac{s_{\mathrm{h}}\left(\mathit{\omega }\right)}{s}}{F}_{\mathrm{s}};\\ & s_{\mathrm{h}}\left(\mathit{\omega }\right)={\displaystyle \frac{s}{\mathrm{1}+{\left(i\mathit{\omega }\mathit{\tau }+{\left(l_{\mathrm{h}}k\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}}},\end{array}\end{array}$$

where h is for “horizontal inhomogeneity” as in Lovejoy et al., 2021; there is an analogous calculation for the FEBE with h≠ $\mathrm{1}/\mathrm{2}$. In the North et al. (1983) 1D model, the top-of-the-atmosphere forcing is exactly a cosine variation, i.e., with a single wavenumber k= 1 cycle around the Earth. The only differences are that the curvature of the Earth was neglected and the Earth's transport properties were assumed to be constant. I nevertheless use Eq. (61) as an approximation for the horizontal transport.

From the data in Table 1 of North et al. (1981), we may deduce the following.

$$\begin{array}{}\text{(62)}& \begin{array}{ll}{F}_{\mathrm{s}}=\left(\mathrm{212}\pm \mathrm{28}\right)e^{-\mathrm{3.27}i}\sin \mathit{\theta };& \mathrm{W}/{\mathrm{m}}^{\mathrm{2}}\\ Q_{\mathrm{s},\mathrm{rad}}=\mathrm{38}{e}^{-\mathrm{3.65}i}\sin \mathit{\theta };& \mathrm{W}/{\mathrm{m}}^{\mathrm{2}}\\ T_{\mathrm{s}}=\mathrm{15.5}{e}^{-\mathrm{3.70}i}\sin \mathit{\theta };& \mathrm{K}\end{array}\end{array}$$

The forcing F_s is the product of the solar constant with the co-albedo (the albedo is 1 minus the co-albedo), θ is the latitude, and the phases are taken with respect to the winter solstice. The variation (about ±13 %) in the amplitude of F_s is due to the latitudinal variation of the co-albedo. In the model, the longwave radiation Q_s,rad and the surface temperature response T_s have exact sinθ dependencies. The phases (in radians) are taken with respect to the winter solstice so that the summer solstice has a phase π= 3.14 rad (in the Northern Hemisphere, 21 June). Due to the co-albedo variations, the actual forcing has a phase of 3.27 rad, peaking on 28 June. Also, the phases of the temperature and longwave emissions are larger at 3.70 rad and 3.65 rad, corresponding to maxima on 26 July and 23 July, respectively (all results are appropriately symmetric for the Southern Hemisphere and for the cold lag following the winter solstice). The nearly identical nature of the phases of temperatures and longwave responses (a 3 d difference, probably not empirically significant) is already support for the model that predicts that they should be in phase. Also note that these lags (of 28, 25 d) are considerably shorter than the 46 d lag (12 August) that would have been obtained had we applied Brunt's heat conductive forcing.

Now use these data to estimate the climate sensitivity, relaxation time τ, and horizontal conduction term l_hk by using the following.

$$\begin{array}{}\text{(63)}& \begin{array}{rl}& s={\displaystyle \frac{T_{\mathrm{s}}}{Q_{\mathrm{s},\mathrm{rad}}}}=\mathrm{0.41}+\mathrm{0.02}i\approx \mathrm{0.41}\mathrm{K}/(\mathrm{W}/{\mathrm{m}}^{\mathrm{2}});\\ & s_{\mathrm{h}}\left(\mathit{\omega }\right)={\displaystyle \frac{T_{\mathrm{s}}}{F_{\mathrm{s}}}}=\left(\mathrm{0.068}\pm \mathrm{0.009}\right)+\left(\mathrm{0.031}\pm \mathrm{0.004}i\right)\mathrm{K}/(\mathrm{W}/{\mathrm{m}}^{\mathrm{2}});\\ & i\mathit{\omega }\mathit{\tau }+{\left(l_{\mathrm{h}}k\right)}^{\mathrm{2}}={\left({\displaystyle \frac{F_{\mathrm{s}}}{Q_{\mathrm{s},\mathrm{rad}}}}-\mathrm{1}\right)}^{\mathrm{2}}=\left(\mathrm{13.20}\pm \mathrm{4.6}\right)\\ & +\left(\mathrm{17.3}\pm \mathrm{5.1}\right)i\end{array}\end{array}$$

From this (with ω= 2π per year), we obtain

$$\begin{array}{}\text{(64)}& \begin{array}{rl}& \mathit{\tau }=\mathrm{2.75}\pm \mathrm{0.8}\mathrm{years}\\ & l_{\mathrm{h}}k=\mathrm{3.63}\pm \mathrm{0.64}.\end{array}\end{array}$$

The relaxation time is within the rough bounds deduced by considering the atmosphere–ocean coupling timescale (≈ 2 years; Hébert et al., 2021), low-frequency climate records (≈ 4.7 ± 2.3 years; Procyk et al., 2020), and high-frequency EBE relaxation times ≈ 4.1 ± 1.1 years (Geoffroy et al., 2013). The ratio of the storage to transfer is $\mathrm{17.3}/\mathrm{13.2}\approx$ 1.3 so that most of the heat is indeed stored and the above homogeneous theory is plausible. The nondimensional l_hk characterizes the typical horizontal divergence over the period of a year. Rather than interpreting it deterministically in terms of a global-scale horizontal variation over a homogeneous Earth, I consider it a nondimensional empirical parameter that will be clarified in future work. In any case, the horizontal transport and storage are in quadrature so that the effect of the transport on the magnitude of sensitivity is smaller, $\left|{\left(i\mathit{\omega }\mathit{\tau }\right)}^{\mathrm{1}/\mathrm{2}}+\mathrm{1}\right|/\left|{\left(i\mathit{\omega }\mathit{\tau }+{\left(l_{\mathrm{h}}k\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}+\mathrm{1}\right|\approx \mathrm{0.88}$ (i.e., about 12 %), but the change in the phase is more substantive (≈ 15 d). Note that the EBE h= 1 value (ignoring transport, with τ= 2.75 years) gives 87 d, i.e., a maximum on 21 September, which is much too late (Fig. 7).

The static climate sensitivity s should be purely real. Its imaginary part is indeed small; it corresponds to 3 d and is probably within the error of the model and empirical estimates, so it will be ignored below. The variable s can be converted to K $/$ (CO₂ doubling) by multiplying it by the canonical value of 3.71 W/m²/(CO₂ doubling) to yield 1.51 K $/$ (CO₂ doubling), which is at the lower end of the IPCC 90 % confidence range (3 ± 1.5 K $/$ (CO₂ doubling)). Since both the methodology and the empirical parameter estimates could be updated and improved, the result is encouraging. In the future, instead of assuming latitudinal constancy with a sinusoidal latitudinal dependence, gridded data could be used and the horizontal conduction approximation (the l_hk term) could be improved.

4 Conclusions

This first paper of two parts proposes a new 2D energy balance equation for macroweather scales: 10 d and longer. It follows the classical energy balance models pioneered by Budyko (1969) and Sellers (1969) and assumes that the dynamics can be adequately modeled by the continuum mechanics heat equation – by advection and diffusion. As reviewed in McGuffie and Henderson-Sellers (2005) and North and Kim (2017), the classical models treat the parts of the atmosphere and ocean that radiatively interact with outer space as a zero-thickness, two-dimensional surface. The complex radiative processes that occur in the vertical direction are only treated implicitly. The dimensionality is then further reduced by zonal averaging.

While this original time-independent model may be reasonable for the long-term (time-invariant) climate states, it is inadequate for treating time-varying anomalies. The key improvement in realism was by made explicitly introducing a vertical coordinate z. Yet, when this was done, it turned out that a detailed vertical model was still unnecessary: all that was required was the existence of a surface layer whose thickness was of the order of the diffusion depth. This is where most of the energy storage occurs, and it determines the vertical temperature derivative at the surface and hence the vertical conductive heat flux. This sensible heat flux is the crucial link between the local radiative imbalances that drive the system, the heat that is stored, and the heat that is transported horizontally. Whereas the Budyko–Sellers models have zero thicknesses, the model has a finite but possibly small thickness; it need only be thick enough to account for energy storage and to determine the surface vertical temperature derivative.

In this first part, I considered only homogeneous zero-dimensional models. These are completely classical, yet as far as I know, they have not been solved with conductive–radiative (linearized) boundary conditions. Using standard Laplace and Fourier techniques, I solved the full depth–time heat equation and showed that its Green's function was identical to a half-order fractional differential equation that directly gives the surface temperature. Although half-order derivatives have occasionally been used in the context of the heat equation (at least since Oldham and Spanier, 1972, 1974; Babenko, 1986), the resulting half-order energy balance equation (the HEBE) is apparently new. Mathematically, the result is a direct consequence of the heat equation, the semi-infinite medium, and conductive–radiative surface boundary conditions.

The consequences are surprisingly far-reaching. For example, the familiar integer-ordered differential equations have exponential Green's functions and short memories. In contrast, the more general fractional-ordered equations such as the HEBE have Green's functions that are generalized exponentials based on power laws and long memories. A general consequence is that while the HEBE respects Newton's law of cooling – i.e., that heat fluxes across a surface are proportional to temperature differences – the relationship between this heat flux and the surface temperature is quite different: it involves a half-order derivative rather than a first-order one. The energy stored is no longer instantaneously determined by the surface temperature, but rather by the entire prior forcing history. Irrespective of the details, we thus expect Earth heat storage processes to generally have long memories.

I also obtained general results for the Earth's response to periodic forcings. Ever since Brunt (1932), Fourier techniques have used the heat equation to model the Earth's temperature response when subjected to a diurnal heat flux forcing. I extended this from the weather regime to the macroweather regime and from diurnally periodic heat forcing to annually periodic radiative–conductive forcing. An immediate consequence is that the surface thermal impedance – equal to the climate sensitivity – is a complex number whose phase determines the lag between the maximum of the forcing (shortly following the summer solstice) and the temperature maximum. Using a simple latitudinally averaged model with empirical parameters, I estimated this complex climate sensitivity and showed how this could readily account for the observed 22–25 d lag, estimating the (static) climate sensitivity at s≈ 0.41 K $/$ (W/m²) with a relaxation time τ≈ 2.75 years.

In Part 2, I extend these zero-dimensional results to the horizontal. I first continue to use Laplace and Fourier techniques to treat the case of homogenous Earth parameters, but with inhomogeneous forcing. Then, with the help of Babenko's method, this is extended to the full inhomogeneous problem with horizontally varying relaxation times, diffusivities, specific heats, climate sensitivities, and forcings.