The Half-order Energy Balance Equation, Part 1: The homogeneous HEBE and long memories

Abstract. The original Budyko–Sellers type 1-D energy balance models (EBMs) consider the Earth system averaged over long times and applies the continuum mechanics heat equation. When these and the more phenomenological zero (horizontal) – dimensional box models are extended to include time varying anomalies, they have a key weakness: neither model explicitly nor realistically treats the surface radiative – conductive surface boundary condition that is necessary for a correct treatment of energy storage. In this first of a two part series, we apply standard Laplace and Fourier techniques to the continuum mechanics heat equation, solving it with the correct radiative – conductive BC's obtaining an equation directly for the surface temperature anomalies in terms of the anomalous forcing. Although classical, this equation is half – not integer – ordered: the Half - ordered Energy Balance Equation (HEBE). A quite general consequence is that although Newton's law of cooling holds, that the heat flux across surfaces is proportional to a half (not first) ordered 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, the relationship is no longer instantaneous. We then consider the case where the Earth is periodically forced. The classical case is diurnal heat forcing; we 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, we show how this can account for the phase lag between the summer maximum forcing and maximum surface temperature Earth response. In part II, we extend all these results to spatially inhomogeneous forcing and to the full horizontally inhomogeneous problem with spatially varying specific heats, diffusivities, advection velocities, climate sensitivities. We consider the consequences for macroweather forecasting and climate projections.



Abstract:
The original Budyko-Sellers type 1-D energy balance models (EBMs) consider the Earth system averaged over long times and applies the continuum mechanics heat equation. When these and the more phenomenological zero (horizontal) -10 dimensional box models are extended to include time varying anomalies, they have a key weakness: neither model explicitly nor realistically treats the surface radiative -conductive surface boundary condition that is necessary for a correct treatment of energy storage.
In this first of a two part series, we apply standard Laplace and Fourier techniques to the continuum mechanics heat equation, solving it with the correct radiative -conductive BC's obtaining an equation directly for the surface temperature anomalies 15 in terms of the anomalous forcing. Although classical, this equation is half -not integer -ordered: the "Half -ordered Energy Balance Equation" (HEBE). A quite general consequence is that although Newton's law of cooling holds, that the heat flux across surfaces is proportional to a half (not first) ordered 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, the relationship is no longer instantaneous. 20 We then consider the case where the Earth is periodically forced. The classical case is diurnal heat forcing; we 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, we show how this can account for the phase lag between the summer maximum forcing and maximum surface temperature Earth response.
In part II, we extend all these results to spatially inhomogeneous forcing and to the full horizontally inhomogeneous problem 25 with spatially varying specific heats, diffusivities, advection velocities, climate sensitivities. We consider the consequences for macroweather forecasting and climate projections.

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 30 possible future climates. The models usually have either zero or one spatial dimension representing respectively the globally or latitudinally averaged meridional temperature distribution (for a review, see e.g. [McGuffie and Henderson-Sellers, 2005 ]).
The fundamental EBM challenge is to model the way that imbalances in incoming short wave and outgoing long wave radiation are transformed into changes in surface temperatures. In an equilibrium climate state, the vertical flux imbalances are 35 transported horizontally. Here we are 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 as well as for surface temperature changes and attendant changes in long wave emissions. EBMs avoid explicit treatment of this critical surface boundary condition, treating it phenomenologically in ways that are flawed; in this two part paper, we show how they can easily be improved with 40 2013], [Sierociuk et al., 2015] and references therein. More generally, fractional derivatives and their equations have a history going back to Leibniz in the 17 th century and their development has exploded in the last decades (for books on the subject, see 75 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 energy budget -has several antecedents arising from the wide range spatial scaling symmetries of atmospheric fields respected by the fluid equations, models and (empirically) by the atmospheric fields themselves (see the reviews [Lovejoy and Schertzer, 80 2013], [Lovejoy, 2019b]). 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], [Hebert, 2017]. However, [Lovejoy, 2019b], [Lovejoy, 2019a], 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 85 fractional energy balance equation (FEBE). Denoting the order of the derivative term in the equation by H, it was shown empirically that if the derivative was of order H ≈ 0.4 -0.5 (rather than the classical EBE value H = 1), that it could account for both the low frequency multidecadal memory [Hebert, 2017] needed for climate projections, as well as the high frequency macroweather (monthly to decadal) memory needed for monthly, seasonal and annual macroweather forecasts, [Lovejoy et al., 2015], [Del Rio Amador and Lovejoy, 2019]. Indeed, the FEBE CRF can be used directly to make climate projections 90 that are compatible with the Coupled Model Intercomparison Project 5 (CMIP5) multi-model ensemble mean projections but with substantially smaller uncertainties (work in progress with R. Procyk). 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 we show 95 how -at least for the H = 1/2 special case-it can be placed on a firmer theoretical basis while simultaneously extending it to two spatial dimensions. Our model is for macroweather temperature anomalies i.e. at time scales longer than the lifetimes of planetary structures, typically 10 days. 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. Our starting point is thus the same as the classical EBMs: radiative, advective and conductive heat transport using the standard continuum mechanics energy 100 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 quite 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, that the surface temperature satisfies the homogeneous HEBE. We relate this to the usual box models, Budyko-Sellers models, and classical diurnal heating models including the notions of thermal admittance and impedance and complex climate 105 sensitivities useful in understanding the annual cycle. We 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 https://doi.org/10.5194/esd-2020-12 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. forcing), and then -using Babenko's method -to the general inhomogeneous case. Part II 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. 110

Conductive and advective heat fluxes
In most of what follows, the earth's spherical geometry plays no role, we use Cartesian coordinates with the z axis pointing upwards and horizontal coordinates x = (x,y) (see however appendix D, part II). We consider that vertical (radiative and conductive), and horizontal (conductive and advective) heat transport occurs on the surface and in the half-volume (x,y,z<0) 115 respectively. Although physically, this means that the atmosphere and ocean are modelled as regions with z≤0, as mentioned, only the vertical surface temperature derivative is ultimately needed and this is well defined if the surface layer is of the order of a few diffusion depths (hundreds of meters). Later, we 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 120 (about 10 days) or longer. Since this is the lifetime of planetary atmospheric structures, much of the actual turbulent atmospheric transport processes are averaged out, giving some justification to the parametrization.
We start with energy transport by diffusion: Fick's law where Qd is the diffusive heat flux vector, k is the thermal diffusivity, r the density, c the specific heat, and T(x,z,t) the temperature. Following standard energy balance models, we use different horizontal ("h") and vertical ("v") coefficients kh(x), kv(x):

125
(1) (the circonflex indicates unit vectors). To include advection, we consider the heat equation for a fluid in a horizontal velocity field vh: Where D/Dt is the advective derivative. The heat equation is therefore: 130 If cr = constant and using the continuity equation, and we can write: (4) Qa is the advective heat flux and T0 is a constant reference temperature (it disappears when the divergence is taken). Taking v = vh, we made the standard assumption that the advection is in the horizontal plane. This is the classical fluid heat equation, it 135 can readily be verified that it conserves energy (integrate both sides over a volume and then use the divergence theorem).
kh(x), kv(x), vh(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.

Radiative heat fluxes 140
At the surface, there is an incoming energy flux : Where F is the anomalous forcing and Q0 (x) is the local solar radiation: S(x) is the local solar constant with local albedo a(x) and the time dependent part of the radiative balance is specified by the 145 additional incoming energy flux, the "forcing" F(x,t). Although in this paper we mostly ignore temporal albedo variations (see however section 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 in place of a(x) in eq. 6 and in place of F in eq. 5. 150 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 (section 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. 155 ∇⋅ cρv is partially balanced by the outgoing that depends on the surface temperature and the effective emissivity e(x): where s is Stefan-Boltzmann constant. The , imbalance drives the system, 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): 160 The derivative term 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 .
If we add the initial condition (or later, ) and the Dirichlet boundary condition at great depth and assume that the system is periodic or infinite in the horizontal, then, in principle, these 165 are enough to determine the temperature for T(x,z<0,t>0) (or eventually, ). Instead of avoiding this conductive -radiative BC below we show how it directly yields an equation for the surface temperature.

The Climatological and anomaly fields
Let us now decompose the heat flux and temperature into time independent (climatological) and time varying (anomaly) components: Qc, Tc and Q, T. As usual, we linearize the outgoing black body radiation, although we do so around the spatially 170 varying surface temperature Tc(x,z = 0) (i.e. not the global average temperature) which yields spatially varying coefficients: (9) (Tc+T is the actual temperature), with climate sensitivity: The linearization is accurate since typical macroweather temperature anomalies are only a few degrees. However due to 175 feedbacks, this formula for the proportionality coefficient -the climate sensitivity -is not accurate; below, we simply consider l(x) to be an empirically determined function of position.
The incoming radiation at the location x drives the system. The radiative imbalance DR going into the subsurface is therefore equal to the conductive flux Qs into the surface; it specifies the conductive-radiative surface boundary condition for Tc and the anomalies T: 180 Where Qd,z is the (upward) vertical component of the heat flux at the surface given by Fick's law: . The conductive -radiative surface boundary conditions for the time independent climate and anomaly temperatures is therefore: (12) 185 l, r, c and k are all presumed to be functions of x. Note: 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, 12 involve the vertical temperature derivative at the surface (i.e. over an infinitesimal layer), lrckv defines the diffusion depth (typically ≈ 10 -100m in thickness, see part II); so that physically the model need only be realistic over this fairly shallow depth where most of the heat is stored. 190 Now, in the temperature eq. 3, replace T by Tc+T. The equation for the time independent climate part is: and for the time-varying anomalies: These equations must now be solved using boundary conditions eqs. 11, 12 for respectively Tc, T and Tc = T = 0 at 195 (all t), and (or see below, ).
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 long wave radiation is approximately linear over the whole range of temperatures (as is commonly assumed in EBMs), this division involves no anomaly smallness assumptions nor assumptions concerning their time averages; the choice of the reference climate depends on the application. Below, we choose anomalies defined in the 200 standard way (although not necessarily with the annual cycle removed, section 3.3), this is adequate for monthly and seasonal forecasts as well as 21 st century climate projections. However, a different choice might be more appropriate for modelling transitions between different climates including possible chaotic behaviours.

The climatological temperature distribution and Budyko-Sellers models 205
In order to simplify the problem, starting with [Budyko, 1969] and [Sellers, 1969], the usual approach to obtaining Tc is somewhat different. First, the climatological temperature field is only defined at z = 0, i.e. Tc(x) = Tc(x,0). Without a vertical coordinate, the climatological radiative imbalance 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). 210
To see how this works, return to eq. 4 for the climatological component and put : (15) (in this formulation, one usually uses the latitude angle instead of the meridional coordinate y see part II, appendix D). The direction of the redirected vertical flux is always away from the equator (y = 0; hence sign(y)), in any event, zonal fluxes will cancel when averaged over latitudinal bands. 215 The usual Budyko-Sellers type models then average Qc over lines of constant latitude yielding a 1-D model:

(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 kh term. There may also be a small advection velocity v but it is not considered to be a true physical velocity but only an ad hoc parameter 220 needed to prevent kh from being negative ( [Sellers, 1969], see also [North et al., 1981]). The horizontal eddy diffusivity kh is often taken as the sum of contributions from water, water vapor and air. In the pure Budyko version, there is no eddy diffusivity, the heat flux is assumed to be proportional to the temperature difference with respect to a reference (e.g. mean) value; . Comparing this with eq. 4 for Qa, we see that this implies that Budyko horizontal heat fluxes are purely advective. 225 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. 230 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 and the corresponding transitions.
Note that the term corresponds to energy storage; in the time dependent case there is no storage.

The nondimensional anomaly equations
No matter how the climate temperature equation is solved, the equation for the time dependent anomaly temperature remains eq. 14. We now rewrite it in a way that brings out the critical mathematical properties. Since rc and kv are only functions of x, eq. 14 can be rewritten: 240 https://doi.org/10.5194/esd-2020-12 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License.
Where we have defined an effective diffusion velocity vd 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, we note that if we include the 245 boundary conditions, the anomaly temperature is entirely determined by the dimensional quantities k, l, r and c. From these, there exists a unique dimensional combination t(x) with dimensions of time, we will see that this controls the relaxation of the system back to thermodynamic equilibrium, it is a "relaxation time". Using kv yields: where lv(x) is the vertical relaxation length of the surface energy balance processes. In part II, table 1and section 3.3 we give 250 some rough parameter estimates. We may also define the horizontal diffusion length lh, speed V and nondimensional advection vector a: For global (zero dimensional) models, t has been estimated as 2 -5 years which is comparable to the classical exponential relaxation time scales mentioned above ( [Hebert, 2017], work in progress with R. Procyk), and in section 3.3 we estimate t ≈ 255 2.75 years.
In order to understand the classical origin of fractional derivatives, it is helpful to consider the homogeneous Seller-type (diffusive transport) heat equation where t, lv and lh are constants and can thus be used to nondimensionalize the operators. t is therefore in units of relaxation times, x in terms of diffusion lengths lh and z in units of diffusion depths lv. By taking l = 1, we effectively have a forcing F with dimensions of temperature. In part I, we consider only the "zero dimensional" equation 260 where the "zero" refers to the number of horizontal dimensions (i.e. only vertical, z and time t). We use the following notation: the first argument is t then horizontal space, then a semicolon followed by the depth z. Circonflexes denote Laplace (time) and Laplace-Fourier (time and horizontal space) transforms.
With z, t in dimensionless form, the homogeneous zero dimensional heat equation is: The transfer is confined to the semi-infinite region z≤0 with boundary conditions: (bottom). At the top (z = 0), the system is forced by the conductive -radiative surface boundary condition: 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 section 3.3 for Fourier methods). 270 Performing a Laplace transform ("L.T.") of the heat equation we obtain: Where the circonflex indicates the Laplace transform in time (with conjugate variable p). Solving: Where A, B are determined by the BC's. Since we require the temperature at depth (z<<0) to remain finite, we must have B = 275 0, hence: To determine A(p), we Laplace transform the surface boundary condition: yielding: 280 (28) It is more convenient to determine the response to the impulse forcing F(t) = d(t); the impulse Green's function.
The above assumes that the subsurface is infinitely deep. If instead it has a finite thickness L, and we take the bottom boundary 285 condition as (rather than ), then and 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 ≈ 100m in the ocean, ≈ 10m for land), the semi-infinite geometry assumption is unimportant. In the following, we therefore ignore finite thickness corrections. 290 Taking the inverse Laplace transform of eq. 29 we obtain the integral representation: (30) (z≤0; where we have used contour integration on the Bromwich integral).

The surface
For the surface, the integral (eq. 30) can be expressed with the help of higher mathematical functions: 295 (31) is the H = 1/2 Mittag-Leffler function (sometimes called a "generalized exponential", also denoted G0,1/2, the "0" for 0 th integral of the impulse, the "1/2" for the order of the derivative for its equation, see below).
For long times after an impulse, the response (t≫1, eq. 33 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) = Q(t) (= 300 1, for t>0, = 0 for t≤0) after which the system eventually attains a new thermodynamic equilibrium. Since , we have the step response (also denoted G1,1/2, eq. 32), and (eq. 33) i.e., a slow power law approach to thermodynamic equilibrium. Figs. 3, 4 show this at different times and depths. With unit step forcing, the boundary condition (eq. 23) indicates that the fraction of the heat flux that is transformed into long wave 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 305 conducted into the subsurface, but that this fraction rapidly vanishes as the surface approaches equilibrium. At equilibrium, the temperature has increased so that the short and long wave fluxes are once again in balance and there is no longer any conductive flux.
For future reference, we give the corresponding step response G1,1/2 = GQ which is the integral of G0,1/2 that describes relaxation to thermodynamic equilibrium when F is a step function. Similarly, the ramp (linear forcing) response G2,1/2 is the integral of 310 the step response, the second integral of the Dirac: For small and large t: The asymptotic equation for the step response (G1,1/2) shows that thermodynamic equilibrium is approached slowly: as t -1/2 . It is this power law step response (with empirical exponent 0.5±0.2) that was discovered semi-empirically by [Hebert, 2017], Lovejoy et al., 2017] and was successfully used for climate projections through to 2100. Similarly, ≈ t 0.4 behaviour was used 320 for macroweather (monthly, seasonal) forecasts close to the short time t 1/2 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: where Dt is the nondimensional time over which (for the TCR) the linear forcing acts. Using t = 4 years, and the standard Dt = 70 325 years for the TCR ramp, we find the plausible ratio TCR/ECS ≈ 0.78.

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 section 3.3. If the surface (z = 0) boundary condition is imposed: 330 (34) then there will be vertical surface gradients that imply that heat is conducted through the surface. To obtain the impulse response Green's function, we take and repeating the Laplace transform approach, we obtain A(p) = 1 (eq. 27 with no derivative term). This yields the following Laplace Transform pairs for the impulse and step Green's function: 335 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 BC's. For example, integrating eq. 30 with respect to time and simplifying, we obtain: 340 Since the step response GQ describes the approach to thermodynamic equilibrium, ( 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. 36 is so that as the figure shows, although they both slowly approach each other and eventually attain thermodynamic equilibrium, that the differences are important (especially in the 345 diffusion layer, z ≈ <1) and they decay very slowly with time and depth, we discuss this further in section 3.3.

Surface temperatures, Fractional derivatives and the HEBE
Let us now introduce the H th order fractional derivative to represent the fractional derivative order H of an arbitrary function f over the domain from t0 to t: (37) 350 Fractional derivatives of order H are most commonly interpreted in the Riemann-Liouville sense ([Podlubny, 1999]) defined by t0 = 0 in the above. This Riemann-Liouville (R-L) formulation is useful for forcings that start at t = 0 (i.e. F = 0 for t≤0), see [Miller and Ross, 1993], [Podlubny, 1999]. The R-L definition is convenient for deterministic forcings, however it singularizes t = 0 whereas we often wish to include periodic or statistically stationary internal stochastic forcings so that (or in the periodic case, the mean over a cycle = 0) is more convenient, in which case we take and 355 hence =0 (or periodic). As discussed in [Lovejoy, 2019a], 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 . The R-L and Weyl fractional derivatives have the following transformation properties: 360 Where w is the Fourier conjugate to t, (see e.g. [Miller and Ross, 1993], [Podlubny, 1999]). In this part I (except for section 3.3), we consider deterministic forcings, putting t0 = 0 in eq. 37, we using (H = 1/2 in eq. 38), 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 ). 370

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", GEBMs (see the reviews [McGuffie and Henderson-Sellers, 2005 ]) whereas in the latter, they may represent the balance across the surface of a homogeneous subsection, a "box". The boxes have spatially uniform 375 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 modelled 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 380 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/area (= S) that is important so that the NLC can be written: S is the heat in the body and Q is the heat flux across the surface into the body (see fig. 6). Teq is the equilibrium temperature, and Z is a transfer coefficient sometimes called 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 Teq (t) = lF(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: 390 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, 42, we may also conclude that thermal impedance Z = l.
While the HEBE and box models obey the NLC, their relationships between the surface heat flux Qs = dS/dt and the surface 395 temperature T are quite different. For example, for forcings starting at time t = t0 , 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., 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 400 comparison, box models satisfy: Where L is the effective thickness of the surface layer and C is the specific heat per area, tbox is the classical EBE relaxation time. [Geoffroy et al., 2013] used a two box model to fit outputs of a dozen GCM and found tbox ≈ 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 2 box model with 405 relaxation scales tbox = 8.4 and 409 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 we take T(t0) = 0 (this is equivalent to fixing the reference of our anomalies), then integrating by parts: 410 (46) Putting H = 1 yields the simple: so that S1 = Sbox.
Over the interval t0 to t, the fractional derivative of order H is defined as the ordinary derivative of the 1-H order fractional integral: (47) 415 Therefore S1/2 = Sbox and:

(48)
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, 2019b], 420 [Lovejoy, 2019a]). 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 = 1/2 not H = 1. Note that when we consider box interfaces with 425 purely conductive heat exchanges (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. dS1/dt = (C/l)dT/dt) can be derived from the heat transport equation. 430 Returning to the nondimensional boundary condition ( ) it is easy to verify, that in order to recover H = 1 heat transfer, one must instead use . We therefore conclude that box model H = 1 transfer is not simultaneously compatible with 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 435 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 1-D and even with this restriction, their time dependent versions have derivatives of the wrong order (see the discussion in part II). 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 reducing it from the phenomenological value H = 1, to H = 1/2. As a bonus, in part II we see that the HEBE can easily be extended from 440 zero to two spatial dimensions, enlarging the scope of energy balance models while simultaneously eliminating these weaknesses.

Conductive versus conductive -radiative boundary conditions
Up until now, we have discussed forcing that is "turned on" at t = 0, this allowed for convenient solutions using Laplace 445 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. We already mentioned that [Brunt, 1932] also considered step function heat forcing, that he claimed might be a plausible model of the diurnal cycle near 450 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 II) 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 p/4 = 3 hours. 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 days which is generally too long (we 455 shall see that it corresponds to an infinite relaxation time).

∂T ∂z z=0
Following [Brunt, 1932] and [Jaeger and Johnson, 1953], let us 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 1-D dimensional heat equation (time and depth, i.e. the dimensional version of eq. 22), we find that the variation of amplitude with depth is: 460 Where Ts is the amplitude of the surface temperature oscillations, it depends on the nature of the forcing, here on the boundary conditions ("s" for "surface"). Following Brunt, using the classical heat surface heat forcing as the surface boundary condition (with this forcing, Fs = Qs is the heat crossing the surface entering the system in the downward direction, see figs. 1, 6) we find: 465 (51) ("heat" for heat forcing), we obtain: Where, Z(w) is the complex frequency dependent thermal impedance, the reciprocal of the thermal admittance. For a given surface heat flux, Z(w) quantifies the surface temperature response (we have written the impedance with the help of l in order 470 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(w), we can thus easily understand the key [Brunt, 1932], [Jaeger and Johnson, 1953] result: that arg(Z(w)) = arg(i -1/2 ) = -p/4 ("arg" indicates the phase).
So far, this approach has only been applied to weather scales (the diurnal cycle). Let's now apply the same approach but with 475 an eye to longer macroweather timescales, notably the annual cycle. The climate sensitivity is an emergent macroweather quantity that is determined by numerous feedbacks that over the weather scales are quite nonlinear but over macroweather scales are considerably averaged (and at least for GCMs, [Hébert and Lovejoy, 2018]) are already fairly linear. 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 yields: 480 Where here Fs is the radiative (downward) forcing radiative flux and Qs and Qs,rad are the surface conductive (into the subsurface) and long wave radiative emission (away from the surface) fluxes respectively. Solving, we obtain the same depth dependence (eq. 50), but with the amplitude of the surface oscillations given by: (54) 485 Where we have introduced the complex climate sensitivity l(w) which by definition is equal to the complex thermal impedance Z(w). In the context of the Earth's energy balance, it is more useful to think in terms of sensitivities than impedances so that below we use l(w). With this, we obtain: Since Arg(i 1/2 ) = p/4 (= 45 o ), we see that as mentioned earlier, the conductive and long wave radiative fluxes are out of phase 490 by 45 o , but the phase of the temperature lags the forcing by Arg(l(w)), which only reaches 45 o in the large t limit (see fig. 7).
Note that we could have deduced eq. 55 directly by Fourier analysis of the HEBE using , but the above allowed us to compare the results with the classical model. The Fourier method allows us to extend the complex climate sensitivity to the more general FEBE: (56) 495 the usual EBE is the H = 1 special case.

Empirical estimates of complex climate sensitivities
Figs. 7, 8 compare the phases and amplitudes of l(w) for the classical and conductive -convective boundary conditions (H = 1/2) HEBE as well as the H = 1 EBE. The plots use w = 2p rad/yr. From fig. 7, we see that taking the empirical value of t in the range 2 -5 years, that the HEBE lag is a little over a month, a result that is close to the observed lag between the summer 500 solstice and maximum temperatures over most land areas. The heat forcing result (p/4 = 1.5 months = 46 days), is already too long for most of the globe and the H = 1 EBE result (close to 3 months = 91 days) 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.
This delays the summer solstice forcing maximum over the ocean, potentially explaining the extra ocean 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 505 by using the latitudinally averaged model discussed in the review [North et al., 1981]. The model uses a 2 nd order Legendre polynomial to take into account the latitudinal variations and a sinusoidal annual cycle with empirically fit parameters that effectively latitudinally average over land and ocean. Empirical parameters are given for the albedo, top of the atmosphere insolation, temperature and outgoing IR radiation.
Before continuing, recall that the zero-dimensional theory discussed here assumes that all radiative flux imbalances are all 510 stored, it ignores horizontal heat transport which -due to the meridional gradients -may be signficant. Although at least for temperature anomalies, we argue that this effect is mostly important at small scales, the magnitude of horizontal transport is not well known and is presumably quite variable from place to place depending on (inhomogeneous) local horizontal transport parameters (see part II). 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 inhomogenous forcing. In part II, we show 515 that for a horizontal wavenumber k, the effect of horizontal transport is to modify the storage term as , therefore for pure periodic horizontal forcing: (57) ("h" for "horizontal inhomogeneity). In North et al's 1-D 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 we neglected the curvature of the 520 Earth and assumed that the Earth's transport properties are constant. We nevertheless use eq. 57 as an approximation for the horizontal transport.
From the data in table 1 of [North et al., 1981] , we may deduce: Where the forcing Fs is the product of the solar constant with the co-albedo (= 1-albedo) and q is the latitude and the phases 525 are taken with respect to the winter solstice. The variation (about ±13%) in the amplitude of Fs is due to the latitudinal variation of the coalbedo. In the model, the long wave radiation Qs,rad and the surface temperature response Ts have exact sinq dependencies. The phases (in radians) are taken with respect to the winter solstice so that the summer solstice has a phase p = 3.14 rads, (in the northern hemisphere, June 21). Due to the coalbedo variations, the actual forcing has a phase = 3.27 rads peaking on June 28th. Also, the phase of the temperature and longwave emissions are larger = 3.70 rad, 3.65 rad corresponding 530 to maxima on July 26 th , July 23 rd respectively (all results are appropriately symmetric for the southern hemisphere and for the cold lag following the winter solstice). The near identity of the phases of temperatures and long wave responses (a three day difference, probably not empirically significant), is already support for the model that predicts that they should be in phase.
We also note that these lags (of 28, 25 days) are considerably shorter than the 46 day lag (Aug 12 th ) that would have been obtained had we applied Brunt's heat conductive forcing. 535 We can use these data to estimate the climate sensitivity, relaxation time t and horizontal conduction term lhk by using the following: From this (with w = 2p/yr), we obtain: Q s,rad = 38e −3.65i sinθ; T s = 15.5e −3.70i sinθ; The relaxation time is within the rough bounds deduced by considering atmosphere -ocean coupling time scale (≈ 2 years, Hebert et al 2020), low frequency climate records (≈ 5 years work with R. Procyk), and the high frequency EBE relaxation times ≈ 4.1±1.1 years [Geoffroy et al., 2013]. We also see that the ratio of the storage to transfer is 17.3/13.2 ≈ 1.3 so that most of the heat is indeed stored so that the above homogeneous theory is plausible. The nondimensional lhk characterizes the typical horizontal transport over the period of a year. Rather than interpreting it deterministically in terms of a global scale 545 horizontal variation over a homogeneous earth, we consider it a nondimensional empirical parameter that we will try to clarify 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: (i.e. about 12%) but the change in the phase is more substantive (≈ 15 days). We can note that the EBE H = 1 value (ignoring transport, with t = 2.75 years) gives 87 days i.e. a maximum on September 21 st which is much too late ( fig. 7). 550 The static climate sensitivity l should be purely real; its imaginary part is indeed small, it corresponds to 3 days and is probably within the error of the model and empirical estimates, it will be ignored below. l can be converted to K/(CO2eq doubling) by multiplying it by the canonical value 3.71 W/m 2 / (CO2eq doubling) to yield 1.51 K/(CO2eq doubling) which is at the lower part of the IPCC 90% confidence range (3±1.5 K/ (CO2eq doubling)). Since both the methodology and the empirical parameter estimates could be updated and improved, the result is encouraging. In future, instead of assuming latitudinal constancy with 555 a sinusoidal latitudinal dependence, gridded data could be used and the horizontal conduction approximation (the lhk term) could be improved.

Conclusions
This first of two parts proposes a new 2D energy balance equation for macroweather scales: ten days and longer. It follows the classical energy balance models pioneered by [Budyko, 1969] and [Sellers, 1969], and assumes that the dynamics can be 560 adequately modelled by the continuum mechanics heat equation -by advection and diffusion. As reviewed in [McGuffie and Henderson-Sellers, 2005 ] 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 565 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 τ = 2.75 ± 0.8yrs 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 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 570 the heat that is transported horizontally. Whereas the Budyko-Sellers models have zero thicknesses, our 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, we considered only homogeneous zero-dimensional models. These are completely classical, yet as far as we know, have not been solved with conductive -(linearized) radiative boundary conditions. Using standard Laplace and Fourier 575 techniques, we solved the full depth-time heat equation and showed that it's 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 [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. 580 The consequences are surprisingly far reaching. For example, the familiar integer ordered differential equations have exponential Green's functions, 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 -that the relationship between this heat flux and the surface temperature is quite different: it involves a half order 585 derivative rather than 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.
We also obtained general results on 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. We 590 extend this from the weather regime to macroweather regime, 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, we estimated this complex climate sensitivity and showed how this could readily account for the observed 22-25 day lag, 595 estimating the (static) climate sensitivity at l ≈ 0.41 K/(W/m 2 ) and relaxation time t ≈ 2.75 years.
In part II, we extend these zero dimensional results to the horizontal. We first continue to use Laplace and Fourier techniques to treat the case of homogenous Earth parameters, but with inhomogeneous forcing. We then -with the help of Babenko's method, extend this to the full inhomogeneous problem with horizontally varying relaxation times, diffusivities, specific heats, climate sensitivities and forcings. 600 https://doi.org/10.5194/esd-2020-12 Preprint. Discussion started: 3 April 2020 c Author(s) 2020. CC BY 4.0 License. Fig. 1: A schematic diagram showing the correct 3D energy balance equations with conductive -radiative surface boundary conditions. Qs is the heat flux across the surface into the subsurface, S is the energy stored in the subsurface per unit surface area.

665
The picture illustrates the thin surface layer (whose thickness is of the order of the diffusion depth, lv with relaxation time t, eq. 20) in which the radiative exchanges between the earth and outer space occur. Teq is the fixed outside temperature, heat will flow as long as the surface temperature Ts ≠ Teq, Z is the thermal impedance (equal here to the climate sensitivity l). To apply the NLC, we need to relate the heat flux to the 690 surface temperature. The lower left shows the consequence of applying heat equation with conductive -radiative BC's, the lower right shows the phenomenological assumption made by box models. The arrows represent heat fluxes, hence the factor l in the denominators. The system is assumed to be horizontally homogeneous and that the subsurface is much thicker than the diffusion depth.

695
Newton's law of cooling  Fig. 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 w = 2p rads/yr) with t in years. The line with short dashes (top) is the usual EBE (H = 1), the solid line is the (H = 1/2) HEBE and the line with long dashes is the classical heat forcing model which is the large t HEBE limit. All curves ignore any net horizontal heat transport. The data analyzed here yield t ≈ 2.75±0.8 years but the actual phase is 700 somewhat shorter due to horizontal heat transport.