the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# 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.

- Article
(1917 KB) - Companion paper
- BibTeX
- EndNote

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 (d*S* $/$ d*t*) is
proportional to the difference between the surface temperature (*T*) and its
long-term equilibrium value (*T*_{eq}), d*S* $/$ d*t* $\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: d*S* $/$ d*t*∝ d*T* $/$ d*t*. 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.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\phantom{\rule{0.125em}{0ex}}\mathit{<}\phantom{\rule{0.125em}{0ex}}\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*
The circumflex indicates unit vectors. To include advection, consider
the heat equation for a fluid in a horizontal velocity field *v*_{h}:

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

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

*Q*_{a} is the advective heat flux and *T*_{0} 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}(

**), and**

*x*

*v*_{h}(

**) are taken to be independent of**

*x**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*_{↓},

where *F* is the anomalous forcing and *Q*_{0}(** x**) is the
local solar radiation:

*Q* is the mean top-of-the-atmosphere flux (≈ 341 W m^{2}), *S*(** x**) is the dimensionless local solar constant with local co-albedo

*a*(

**) (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”**

*x**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*(

**) 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**

*x**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**):

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*):

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.

## 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:

where *T*_{c}+*T* is the actual temperature. With climate sensitivity, it yields

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*:

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(\right)}_{\frac{\partial T}{\partial z}}z=\mathrm{0}$. The conductive–radiative
surface boundary conditions for the time-independent climate and anomaly
temperatures is therefore

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

And the equation for the time-varying anomalies is

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}(

**, 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).**

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

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:

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:

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

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

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

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

**:**

*α*
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

where *ζ* is the dimensionless horizontal transport operator. The reference temperature *T*_{0} 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

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:

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^{6}J $/$ (m^{3}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_{2}doubling value 3 ± 1.5 K, a 90 % confidence interval, and 3.71 W m^{2}for CO_{2}doubling, the global mean value is*s*≈ 0.8 ± 0.4 K W m^{2}, with regional values a factor of ≈ 2 higher or lower (IPCC AR5), yielding*ρ**c**s*≈ 3 × 10^{6}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^{8}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^{6}Wm^{−1}/K from which we estimate the horizontal (eddy) diffusivity as*κ*_{h}=*K*/(*ρ*c) ≈ 1 m^{2}/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^{−4}m^{2}/s are typical (Houghton et al., 2001). For soil, rough values are*κ*_{v}≈ 10^{−6}m^{2}/s (wet) and*κ*_{v}≈ 10^{−7}m^{2}/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^{8}s to obtain*κ*_{v}≈ 10^{−5}m^{2}/s, which is close to the model values. - f.
Diffusion depth

*l*_{v}. Using*l*_{v}=*κ*_{v}*ρ**c**s*, we find for the ocean and soils*l*_{v}≈ 300 m and ≈ 3–10 m, respectively. Using the global estimates,*κ*_{v}≈ 10^{−5}–10^{−4}m^{2}/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^{2}/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}–3 × 10^{−4}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^{2}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^{−5}–10^{−4}m/s (Eq. 4); hence, using*V*≈ 10^{−4}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

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

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

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

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

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

yielding

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

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:

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.

*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.

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.

For small and large *t*, we have the following.

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,

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.

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.

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.

### 3.1.6 Surface temperatures, fractional derivatives, and the HEBE

Let us now introduce the *h*th-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*_{0} to *t*:

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}= 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:

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}= 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.

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:

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

*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^{2} K/W); its
reciprocal *Y* is the surface “thermal admittance” (see the next section).
Identifying the equilibrium temperature with *T*_{eq}(*t*)= *s**F*(*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:

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*.

While the HEBE and box models both obey the NLC, the relationships between
the surface heat flux *Q*_{s}= d*S* $/$ d*t* and the surface temperature *T* are quite
different. For example, for forcings starting at time *t*=*t*_{0 }, using the HEBE we have

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

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 *h*th-order heat storage as

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

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

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

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

Combining this with the NLC, in both cases we obtain

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$)d*T* $/$ d*t*) can be derived from the heat transport
equation. Returning to the nondimensional boundary condition $\left({\left(\right)}_{\frac{\partial T}{\partial z}}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

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*_{s}*e*^{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

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

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*_{s}*e*^{iωt} yields

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

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

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:

The usual EBE is the *h*= 1 special case.

### 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,

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.

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*_{h}*k* by using the following.

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

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*_{h}*k* 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_{2} doubling) by multiplying it by the canonical value of 3.71 W/m^{2}/(CO_{2} doubling) to yield 1.51 K $/$ (CO_{2} doubling), which is at the lower end of
the IPCC 90 % confidence range (3 ± 1.5 K $/$ (CO_{2} 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*_{h}*k* term)
could be improved.

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^{2}) 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.

The figures were produced using a standard numerical integration package.

The data in Sect. 3.3.2 were taken from the cited source.

The author declares that there is no conflict of interest.

I acknowledge discussions with Lenin Del Rio Amador, R. Procyk, Raphael Hébert, Dave Clarke, and Cécile Penland.

This paper was edited by Anders Levermann and reviewed by Peter Ashwin and one anonymous referee.

Atanackovic, M., Pilipovic, S., Stankovic, B., and Zorica, D.: Fractional Calculus with applications in mechanics: variations and diffusion processes, Wiley, London, UK, 313 pp., 2014.

Babenko, Y. I.: Heat and Mass Transfer, Khimiya, Leningrad, 1986 (in Russian).

Baleanu, D., Diethelm, K., Scalas, E., and Trujillo, J. J.: Fractional Calculus: Models and Numerical Methods,, World Scientific, Singapore, 400 pp., 2012.

Blender, R. and Fraedrich, K.: Long time memory in global warming simulations, Geophys. Res. Lett., 30, 1769–1773, 2003.

Brunt, D.: Notes on radiation in the atmosphere, Q. J. Roy. Meterol. Soc., 58, 389–420, 1932.

Budyko, M. I.: The effect of solar radiation variations on the climate of the earth, Tellus, 21, 611–619, 1969.

Bunde, A., Eichner, J. F., Kantelhardt, J. W., and Havlin, S.: Long-term memory: a natural mechanism for the clustering of extreme events and anomalous residual times in climate records, Phys. Rev. Lett., 94, 1–4, https://doi.org/10.1103/PhysRevLett.94.048701, 2005.

Byrne, G. F. and Davis, J. R.: Thermal inertia, thermal admittance, and the effect of layers, Remote Sens. Environ., 9, 295–300, https://doi.org/10.1016/0034-4257(80)90035-8, 1980.

Del Rio Amador, L. and Lovejoy, S.: Predicting the global temperature with the Stochastic Seasonal to Interannual Prediction System (StocSIPS), Clim. Dynam., 53, 4373–4411, https://doi.org/10.1007/s00382-019-04791-4, 2019.

Del Rio Amador, L. and Lovejoy, S.: Long-range Forecasting as a Past Value Problem: Using Scaling to Untangle Correlations and Causality, Geophys. Res. Lett., https://doi.org/10.1029/2020GL092147, 2021a.

Del Rio Amador, L. and Lovejoy, S.: Using regional scaling for temperature forecasts with the Stochastic Seasonal to Interannual Prediction System (StocSIPS), Clim. Dynam., 1432–0894, https://doi.org/10.1007/s00382-021-05737-5, 2021b.

Donohoe, A., Dawson, E., Mcmurdie, L., Battisti, D. S., and Rhines, A.: Seasonal asymmetries in the lag between insolation and surface temperature, J. Climate, 33, 3921–3945 https://doi.org/10.1175/jcli-d-19-0329.1, 2020.

Dwyers, H. A. and Petersen, T.: Time-dependent energy modelling, J. Appl. Meteorol., 12, 36–42, 1975.

Geoffroy, O., Saint-Martin, D., Olivié, D. J., Voldoire, A., Bellon, G., and Tytéca, S.: Transient climate response in a two-layer energy-balance model. part i: Analytical solution and parameter calibration using cmip5 aogcm experiments, J. Climate, 26, 1841–1857, 2013.

Ghil, M. and Lucarini, V.: The physics of climate variability and climate change, Rev. Mod. Phys., 92, 035002, https://doi.org/10.1103/RevModPhys.92.035002, 2020.

Gottwald, G. A., Crommelin, D. T., and Franzke, C. L. E., Stochastic Climate Theory, in: Nonlinear and Stochastic Climate Dynamics, edited by: Franzke, C. L. E. and O'Kane, T. J., Cambridge University Press, 209–240, https://doi.org/10.1017/9781316339251.009, 2017.

Hahn, D. W. and Ozisk, M. N.: Heat Conduction, 3rd edition ed., Wiley, New York, 2012.

Hebert, R.: A Scaling Model for the Forced Climate Variability in the Anthropocene, MS thesis, McGill University, Montreal, 2017.

Hébert, R. and Lovejoy, S.: The runaway Green's function effect: Interactive comment on “Global warming projections derived from an observation-based minimal model” by K. Rypdal, Earth Syst. Dynam. Disc., 6, C944–C953, 2015.

Hébert, R. and Lovejoy, S.: Regional Climate Sensitivity and Historical Based Projections to 2100, Geophys. Res. Lett., 45, 4248–4254, https://doi.org/10.1002/2017GL076649, 2018.

Hébert, R., Lovejoy, S., and Tremblay, B.: An Observation-based Scaling Model for Climate Sensitivity Estimates and Global Projections to 2100, Clim. Dynam., 56, 1105–1129, https://doi.org/10.1007/s00382-020-05521-x, 2021.

Hilfer, R.: Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000.

Houghton, J. T., Ding, Y., Griggs, D. J., Noguer, M., van der Linden, P. J., Dai, X., Maskell, K., and Johnson, C. A.: Climate Change 2001: The Scientific Basis, Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, 2001.

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, 2013.

Jaeger, J. C. and Johnson, C. H.: Note on diurnal temperature variation, Pure Appl. Geophys., 24, 104–106, 1953.

Klafter, J., Lim, S., and Metzler, R.: Fractional Dynamics: Recent Advances, World Scientific, Singapore, 2012.

Lovejoy, S.: Using scaling for macroweather forecasting including the pause, Geophys. Res. Lett., 42, 7148–7155 https://doi.org/10.1002/2015GL065665, 2015.

Lovejoy, S.: Weather, Macroweather and Climate: our random yet predictable atmosphere, Oxford University Press, New York, 334 pp., 2019a.

Lovejoy, S.: Fractional relaxation noises, motions and the fractional energy balance equation, Nonlin. Processes Geophys. Discuss. [preprint], https://doi.org/10.5194/npg-2019-39, in review, 2019b.

Lovejoy, S.: The half-order energy balance equation – Part 2: The inhomogeneous HEBE and 2D energy balance models, Earth Syst. Dynam., 12, 489–511, https://doi.org/10.5194/esd-12-489-2021, 2021.

Lovejoy, S. and Schertzer, D.: The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, Cambridge, 496 pp., 2013.

Lovejoy, S., del Rio Amador, L., and Hébert, R.: The ScaLIng Macroweather Model (SLIMM): using scaling to forecast global-scale macroweather from months to decades, Earth Syst. Dynam., 6, 637–658, https://doi.org/10.5194/esd-6-637-2015, 2015.

Lovejoy, S., Del Rio Amador, L., and Hébert, R.: Harnessing butterflies: theory and practice of the Stochastic Seasonal to Interannual Prediction System (StocSIPS), in: Nonlinear Advances in Geosciences, edited by: Tsonis, A. A., Springer Nature, Switzerland, 305–355, 2017.

Lovejoy, S., Procyk, R., Hébert, R., and del Rio Amador, L.: The Fractional Energy Balance Equation, Q. J. Roy. Meteorol. Soc., https://doi.org/10.1002/qj.4005, 1–25, 2021.

Márquez, J. M. A., Bohórquez, M. A. M., and Melgar, S. G.: Ground Thermal Diffusivity Calculation by Direct Soil Temperature Measurement. Application to very Low Enthalpy Geothermal Energy Systems, Sensors-Basel, 16, 306, https://doi.org/10.3390/s16030306, 2016.

McGuffie, K. and Henderson-Sellers, A.: A Climate Modelling Primer, Third Edition ed., John Wiley & Sons Ltd, Chichester, England, 2005.

Meyer, R. F.: A heat-flux-meter for use with thin film surface thermometers: a report, Report, National Research Council of Canada, Ottawa, 1960.

Miller, K. S. and Ross, B.: An introduction to the fractional calculus and fractional differential equations, John Wiley and Sons, New York, 366 pp., 1993.

Mori, H.: Transport, Collective Motion, and Brownian Motion, Prog. Theor. Phys., 33, 423–455, https://doi.org/10.1143/PTP.33.423, 1965.

Myrvoll-Nilsen, E., Sørbye, S. H., Fredriksen, H.-B., Rue, H., and Rypdal, M.: Statistical estimation of global surface temperature response to forcing under the assumption of temporal scaling, Earth Syst. Dynam., 11, 329–345, https://doi.org/10.5194/esd-11-329-2020, 2020.

North, G. R. and Coakley, J. A.: Differences between seasonal and mean annual energy balance model calculations of climate and climate sensitivity, J. Atmos. Sci., 36, 1189–1204, 1979.

North, G. R., Cahalan, R. F., and Coakley, J. J. A.: Energy balance climate models, Rev. Geophys. Space Phy., 19, 91–121, 1981.

North, G. R., Mengel, J. G., and Short, D. A.: Simple Energy Balance Model Resolving the Seasons and the Continents Application to the Astronomical Theory of the Ice Ages, J. Geophys. Res., 88, 6576–6586, 1983.

North, R. G. and Kim, K. Y.: Energy Balance Climate Models, Wiley-VCH, Weinheim, Germany, 369 pp., 2017.

Oldham, K. B.: Diffusive transport to planar, cylindrical and spherical electrodes, J. Electroanal. Chem. Interfacial Electrochem., 41, 351–358, 1973.

Oldham, K. B. and Spanier, J.: A general solution of the diffusion equation for semi infinite geometries, J. Math. Anal. Appl., 39, 665–669, 1972.

Oldham, K. B. and Spanier, J.: The Fractional Calculus, Academic Press, reprinted by Dover 2006, New York, USA, 1974.

Podlubny, I.: Fractional Differential Equations, Academic Press, San Diego, United States, 340 pp., 1999.

Procyk, R., Lovejoy, S., and Hébert, R.: The Fractional Energy Balance Equation for Climate projections through 2100, Earth Syst. Dynam. Discuss. [preprint], https://doi.org/10.5194/esd-2020-48, in review, 2020.

Rybski, D., Bunde, A., Havlin, S., and von Storch, H.: Long-term persistance in climate and the detection problem, Geophys. Resear. Lett., 33, L06718, https://doi.org/10.1029/2005GL025591, 2006.

Rypdal, K.: Global temperature response to radiative forcing: Solar cycle versus volcanic eruptions, J. Geophys. Res., 117, D06115, https://doi.org/10.1029/2011JD017283, 2012.

Rypdal, K.: Global warming projections derived from an observation-based minimal model, Earth Syst. Dynam., 7, 51–70, https://doi.org/10.5194/esd-7-51-2016, 2016.

Rypdal, K., Rypdal, M., and Fredriksen, H., Spatiotemporal Long-Range Persistence in Earth's Temperature Field: Analysis of Stochastic-Diffusive Energy Balance Models, J. Climate, 28, 8379–8395, https://doi.org/10.1175/JCLI-D-15-0183.1, 2015.

Rypdal, M. and Rypdal, K.: Long-memory effects in linear response models of Earth's temperature and implications for future global warming, J. Climate, 27, 5240–5258, https://doi.org/10.1175/JCLI-D-13-00296.1, 2014.

Sellers, W. D.: A global climate model based on the energy balance of the earth-atmosphere system, J. Appl. Meteorol., 8, 392–400, 1969.

Sierociuk, D., Dzielinski, A., Sarwas, G., Petras, I., Podlubny, I., and Skovranek, T.: Modelling heat transfer in heterogeneous media using fractional calculus, Philos. T. R. Soc. A, 371, 20120146 https://doi.org/10.1098/rsta.2012.0146, 2013.

Sierociuk, D., Skovranek, T., Macias, M., Podlubny, I., Petras, I., Dzielinski, A., and Ziubinski, P.: Diffusion process modeling by using fractional-order models, Appl. Math. Comput., 257, 2–11 https://doi.org/10.1016/j.amc.2014.11.028, 2015.

Stubenrauch, C. J., Chédin, A., Rädel, G., Scott, N. A., and Serrar, S.: Cloud Properties and their seasonal and diurnal variability from TOVS Path-B, J. Climate, 19, 5531–5553 https://doi.org/10.1175/JCLI3929.1, 2006.

Tarasov, V. E.: Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media, Higher Education Press, Beijing, China, 2010.

Trenberth, K. E., Fasullo, J. T., and Kiehl, J.: Earth's global energy budget, B. Am. Meteorol. Soc., 90, 311–323 https://doi.org/10.1175/2008BAMS2634.1, 2009.

van Hateren, J. H.: A fractal climate response function can simulate global average temperature trends of the modern era and the past millennium, Clim. Dynam., 40, 2651, https://doi.org/10.1007/s00382-012-1375-3, 2013.

Varotsos, C. A., Efstathiou, M. N., and Cracknell, A. P.: On the scaling effect in global surface air temperature anomalies, Atmos. Chem. Phys., 13, 5243–5253, https://doi.org/10.5194/acp-13-5243-2013, 2013.

Warren, S. G. and Schneider, S. H.: Seasonal simu-lation as a test for uncertainties in the parameterizations of a Budyko-Sellers zonal climate model, J. Atmos. Sci., 36, 1377–1391, 1979.

West, B. J., Bologna, M., and Grigolini, P.: Physics of Fractal Operators, Springer, New York, USA, 354 pp., 2003.

Zhuang, K., North, G. R., and Stevens, M. J.: A NetCDF version of the two-dimensional energy balance model based on the full multigrid algorithm, SoftwareX, 6, 198–202, https://doi.org/10.1016/j.softx.2017.07.003, 2017.

Ziegler, E. and Rehfeld, K.: TransEBM v. 1.0: Description, tuning, and validation of a transient model of the Earth’s energy balance in two dimensions, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2020-237, in review, 2020.

Zwanzig, R.: Nonlinear generalized Langevin equations, J. Stat. Phys., 9, 215–220, https://doi.org/10.1007/BF01008729, 1973.

Zwanzig, R.: Nonequilibrium Statistical Mechanics, Oxford University Press, USA, 2001.