The Half-order Energy Balance Equation, Part 2: The inhomogeneous HEBE and 2D energy balance models

: In part I, we considered the zero-dimensional heat equation showing quite generally that conductive


Introduction
In part I, we showed that when the surface of a body exchanges heat both conductively and radiatively, that its flux depends on the half order derivative of the surface temperature.This implies that energy stored in the subsurface effectively has a huge power law memory.This contrasts with the usual phenomenological assumption used notably in box models (including zero dimensional global energy balance models) that the order of derivative is an integer (one) and that on the contrary, the memory is only exponential (short).The result followed directly by assuming that the continuum mechanics heat equation was obeyed and the depth of the media was of the order of a few diffusion depths; for the Earth, perhaps several hundred meters.The basic result was a classical application of the heat equation barely going beyond results that [Brunt, 1932] already found "in any textbook".
A consequence was that although Newton's law of cooling is obeyed, that the temperature obeyed the half-order energy balance equation (HEBE) rather than the phenomenological first order Energy balance Equation (EBE).When applied to the Earth, the HEBE and its implied long memory explains the success of both climate projections through to 2100 [Hebert, 2017], [Lovejoy et al., 2017], [Hébert et al., 2020] and macroweather (monthly, seasonal) temperature forecasts [Lovejoy et al., 2015], [Del Rio Amador and Lovejoy, 2019], [Del Rio Amador and Lovejoy, 2020a;Del Rio Amador and Lovejoy, 2020b].
We also considered the responses to periodic forcings showing that surface heat fluxes and temperatures are related by a complex thermal impedance (Z(w), w is the frequency).In the Earth system, Z(w) = s(w) where s(w) is the complex climate sensitivity that we estimated from a simple semi-empirical model.
Although in part 1 we discussed the classical 1-D application of the heat equation to the Earth's latitudinal energy balance (Budyko-Sellers models) -especially their ad hoc treatment of the surface boundary condition -we restricted the discussion to zero horizontal dimensions.In this part II, we first (section 2) extend the part I treatment to systems with homogeneous properties but with inhomogeneous forcings, first in the horizontal plane (section 2.1, 2.2), then -following Budyko-Sellerslatitudinally varying on the sphere (section 2.3).The homogeneous case is quite classical and can be treated with standard Laplace and Fourier techniques, it leads to the (horizontally) Generalized HEBE: the GHEBE.Although the GHEBE has a more complex (space-time) fractional derivative operator that is unlike anything we know of in the literature -like the HEBE, it can nevertheless be given precise meaning via its Green's function.
In section 3, we derive the inhomogeneous GHEBE and HEBE needed for applications.This is done by using of Babenko's method [Babenko, 1986] which is essentially a generalization of Laplace and Fourier transform techniques.The challenge with Babenko's method is to interpret the inhomogeneous space-time fractional operators.Following Babenko, we do this using both high and low frequency expansions corresponding respectively to processes dominated by storage and by horizontal heat transport.The long time limit describes the new energy balance climate state that results when the forcing is increased everywhere and held fixed: for the model this corresponds to equilibrium.We also include several appendices focused on empirical parameter estimates (appendix A), the implications for two point and space-time temperature statistics (when the system is stochastically forced, internal variability, appendix B), and finally (appendix C), the changes needed to account for the Earth's spherical geometry, including the definition of fractional operators on the sphere.

The homogeneous GHEBE
In part I we recalled the heat equation for the time-varying temperature anomalies (T) with diffusive and (horizontal) effective advective velocity (v): (1) (This is written in the still general form of eq. 19, part I).kh, kv are horizontal and vertical thermal diffusivities, z the vertical coordinate (pointing upwards, the Earth is z≤0), t the time, x = (x,y) the horizontal coordinates, (the circonflexes indicate unit vectors).These equations must now be solved using the conductive-radiative surface boundary condition: (2) r, c are the fluid densities and specific heats s is the climate sensitivity and F is the anomaly forcing.The initial conditions are T = 0 at (all t), and (Riemann-Liouville) or below, (Weyl).
In part I, we nondimensionalized the zero-dimensional homogeneous operators by nondimensionalizing time by the relaxation time: (with ) and nondimensionalizing the vertical distance by the vertical diffusion depth: , with .Considering now the full equation with advective and diffusive transport, we nondimensionalize the horizontal coordinates by the horizontal diffusion length: , (with ) and use the nondimensional advection velocity (with speed ).If we now take s = 1 (equivalent to using dimensions of temperature for the forcing F), we obtain: (3) For the heat equation and the conductive-radiative surface boundary condition respectively.For initial conditions such that T = 0 for t≤0, as in part I, we take Laplace transforms in time, but we now take Fourier transforms in the horizontal: (4) Where "F.T." is the Fourier transform in horizontal space, k for the conjugate of x, (the vector modulus) with conjugate variable (as usual, ).Fourier transforms in space are convenient for either infinite horizontal media, or media with periodic horizontal boundary conditions.In appendix C, we consider the changes needed to account for spherical geometry.

When
, the solution and where Gd is the impulse (Dirac) response Green's function, part I, eq.30.From eq. 4, we see that this is the same as the zero dimensional equation (eq.24, part I) but with i.e. for the corresponding Green's function: A note on notation: the first argument is time, with the vertical separated by a semi-colon.When there is a horizontal coordinate it comes after time, before the semicolon.With this notation, the right hand side of eq. 5 is the L.T. of the zero-dimensional (time-depth) Green's function , the left hand side is the Laplace (time) and Fourier transform (horizontal, space) transform.
We can now use the basic Laplace shift property: To conclude that: Decomposing this into a circularly symmetric diffusion part and a factor that shifts phases, we obtain: By circular symmetry of , its inverse (2-D) Fourier transform reduces to an inverse Hankel transform ("H.T."). Using: (9) We therefore obtain for the diffusive part of the surface impulse response (i.e. the response with source spatial forcing ): (10) Where is the zero-dimensional impulse response.If needed, its integral representation is given in eq.34, part I.The last step is to take into account the advective term associated with the phase shift .For this final step, we use the Fourier shift theorem to obtain: This is the general surface result for the diffusive-advective transport part of the spatially homogeneous case.As expected, the advective transport simply displaces the centre of the impulse response with nondimensional velocity a.As usual, the solutions for arbitrary forcing F(t,x) can be obtained by convolution.
For the surface we obtain the simpler expressions: (12) (see eq. 35, part I).From these, the general surface results including advection are obtained with , i.e.
. Since the advection term has this simple consequence, below we take a = 0, considering only diffusive transport, advection can easily be included if needed (i.e.below, we take ).
To better understand the impulse response, fig. 1 shows this surface for various radial distances r and fig. 2 shows the corresponding time dependence of the time integral of Gd; the unit step response GQ for various distances r, illustrating the power law approach to equilibrium at large t (discussed in section 2.2).The corresponding long time, short distance expansions are: Where is the Green's function for the (spatial Dirac) "hotspot" equilibrium response discussed below (eq. 20).Note that the leading term in is independent of r, and the leading term in the approach to equilibrium is also independent of r.
Just as we derived the zero-dimensional HEBE by showing that it had the same Green's function as the z = 0 transport equation Green's function, we can likewise derive the homogeneous Generalized Half-Order Energy Balance Equation (GHEBE) which is the space-time surface equation whose Green's function is given in eq.12.Following the derivation of the HEBE, in part I eq.29, and replacing we obtain: Hence, for z = 0: The left hand equation is the homogeneous GHEBE whose Green's function is given by eq.12.We have therefore found a surprisingly simple explicit formula for the (inverse) half-order space-time GHEBE operator: where " " indicates convolution.This allows us to give a precise interpretation of the half-order operator.Therefore the dimensional, homogeneous, GHEBE and its full solution are: (17) G therm,δ r,0 ( ) ("surf" is the surface over which the forcing acts, the bottom line uses the explicit eq. 12 for Gd).
The above shows that even with the purely classical integer-ordered Budyko-Sellers type heat equation, that surface temperatures already obey long memory, half order equations.However, it is not certain that the classical heat equation is in fact the most appropriate model.Straightforward generalizations to fractional heat equations -where lead directly to fractional energy balance equations for surface temperatures, we investigate fractional heat equations elsewhere.Physically, this generalization from the classical fractional value H = 1/2 could be a consequence of turbulent diffusive transport which since at least Richardson been known to have anomalous diffusion.

Energy balance, equilibrium
If F(t,x) = 0 then there is a radiative energy balance at time t, point x, but the temperature may be changing.However, if F(t,x) = 0 for a long enough time, and for all x, then the time derivatives ( ) vanish and Earth is in a steady energy balance ("climate") state, Tclim(x), so that the temperature anomaly T(t,x) = 0. Now consider a step function increase .Then as , the time derivatives will vanish and a new (steady) climate state (with temperature ) will be reached in which the horizontal transport and anomalous black body emission balance the new forcing: . The new state is steady in time and is in energy balance with outer space and its local surroundings, but it is not strictly correct to describe as one of thermal equilibrium.This is because thermal equilibrium would imply that the temperature everywhere is constant (thermodynamic equilibrium is an even more stringent condition).Nevertheless the term "radiative equilibrium" is commonly used in the context of planetary energy balance, so we will use the terms energy balance and equilibrium synonymously.
Let us now investigate the equilibrium state.Since d/dt = 0, the conjugate variable p = 0, with this and a = 0 in eq.15, we obtain the equation for the (spatial) surface impulse response for equilibrium (subscript "eq"): (18) i.e. the same as eq. 4 but with p = 0 (and a = 0) hence: The equilibrium surface temperature (spatial) impulse (Dirac "hotspot") Green's function is therefore: Where H0 is the zeroth order Struve function and Y0 is the zeroth order Bessel function of the second kind.For large r, we have the expansions: ; r ≈ 0 The 1/r 3 asymptotic decay is fast and implies that spatial hotspots remain fairly localized; indeed, it is easy to show that if instead we had a Dirac surface heat flux source driving the system (i.e. with surface BC i.e. without radiation) that the decay would be the much faster (1/r).Forcing inhomogeneities thus remain much more localized than would otherwise be the case.
To study the convergence to equilibrium, consider a simple model of a surface "hot spot" where the forcing is confined to a unit circle and turned on and held at a constant unit temperature at t = 0.This is the spatial equivalent of a step forcing in space, we combine it with a step (Heaviside) in time: P1(r) is the corresponding indicator function.We now use the transform pair to perform the convolution: (23) (J1 is the first order Bessel function of first kind).Taking the limit we obtain the equilibrium temperature distribution.
Alternatively we could find it directly by from eq. 19: G eq,δ r,0 G eq,δ r;0 Fig. 4 shows the cross section as a function of the distance from the circle's center at various times (the inverse Hankel transforms were done numerically).We note that the temperature rises very quickly at first, then slowly reaches equilibrium (thick).The figure also shows (dashed) the equilibrium when the forcing is purely due to unit conductive heating over the unit circle.The difference between the dashed and the thick equilibrium curves are purely due to the radiative loses in the latter.(Note that in the zero-dimensional case (part I), using pure heating forcing boundary conditions leads to diverging temperatures, there is no equilibrium.This explains why Brunt instead used temperature forcing boundary conditions.Here, in two horizontal dimensions, boundary conditions that impose a fixed temperature over the circle are problematic since they imply infinite horizontal temperature gradients and infinite horizontal heat fluxes).
Figs. 5, 6 shows the same evolution but with temperature as a function of time for various distances (fig.5) and as contours in space-time (fig.6).We see that equilibrium is largely established in the first two relaxation times (here t = 1) and most of the perturbation is confined to two horizontal diffusion distances (here, lh = 1).

Comparison of the HEBE with the standard 1-D Budyko -Sellers model on a sphere
It is helpful to clearly understand the similarities and differences between the HEBE and the usual 1-D (latitudinal) B-S approach (see the comprehensive monograph [North and Kim, 2017], and see [Zhuang et al., 2017], [Ziegler and Rehfeld, 2020] for recent applications, development).Since the latter model is on a sphere but with only latitudinal dependence, we write the horizontal transport term using gradient and divergence operators: with q = colatitude and µ = cos q.In standard notation [North and Kim, 2017]) the B-S equation is thus written: (25 Where C is the specific heat per area, the thermal conductivity per radian of arc, B is the climate feedback parameter, the inverse of climate sensitivity (B =1/s), Q0, the solar constant, H the heat function, S is the insolation distribution function, a is the co-albedo and A is the constant term from the linearization of the black-body emission.If we measure temperatures with respect to the mean (reference) Earth temperature so that the A term balances the mean forcing), then the B-S equation with dimensionless operators can be written: where F is the anomaly with respect to the global average.
In part I section 3.1.1,we expressed the horizontal transport operator in terms of the transport coefficient DF that allows us to write the HEBE in the form: where .Using for the transport operator, we obtain the 1-D HEBE on the sphere: (28) In the case of constant thermal diffusion coefficients we may solve both equations using Legendre polynomials Pn(µ) that are eigenfunctions of the Laplacian: (with boundary conditions at the poles being zero horizontal heat flux, see also appendix C for more general results on the sphere).Expanding the temperature and forcing in terms of the Legendre polynomials and taking Laplace transforms of the coefficients in time, we obtain: (29) We then obtain equations for the Laplace transform of the n th Legendre coefficients: So that: (31) In real space: 225 (32) (note that the generalization to the FEBE is obtained by the replacement so that whereas so that is not a special case of the FEBE).Using: (33) 230 (eq.35, part I), combining this with eq.32, we obtain for the impulse responses: (34) Integrating these with respect to t, we obtain the step responses: ) The long time limit represents Earth energy balance (equilibrium): (36) If x<0, then there is an unphysical divergence so that sDF must be >0.Since Pn(µ) has n zeroes, n plays the role of wavenumber, it specifies structures of horizontal size ≈ pR/n.Therefore we see that the B-S model (where G falls off as n -2 ) will yield a much smoother equilibrium temperature than the HEBE where it falls off as n -1 .Note that when generalized from the HEBE to the FEBE (with p→p 2H ), this equilibrium result is unchanged.
For the HEBE, the short and long time behaviours are: (37) The asymptotic response for is interesting because it tells us how quickly equilibrium is reached.When n = 0 we have P0(µ) = 1, so that this component corresponds to the mean.Since we see that it is identical to the zerodimensional result in part I: equilibrium is approach in a power law fashion (t -1/2 for large t), whereas for n = 0, the B-S model ξ F ,0 = 0 approach to equilibrium is exponential.However for n≥1, HEBE power law terms are exponentially damped, with exponential decay time: whereas the B-S model is exponentially damped for all n with .
In order to make a more detailed comparison between the models, we can follow [North and Kim, 2017] who consider a model with constant DS-B and that is north-south symmetric so that the odd numbered polynomials vanish.They empirically give the climate equilibrium values for n = 0, 2, 4; the (constant) n = 0 term is used to obtain the mean temperature 288K.Other pertinent empirical data are s = 1/B = 0.50 KW -1 m 2 , F2 = -180.7 W/m 2 , F4 = 20.8K,T2 = -30K, T4 = -4K.From eq. 36 for the equilibrium temperature Green's function, we obtain: . The n = 2 relationship is use to estimate = 0.67 Wm -2 K -1 , with this estimate, we obtain ≈ 1.35K which is not far from the empirical estimate T4 = -4K ( [North and Kim, 2017]), it also yields the dimensionless quantity sDB-S = 0.33.If we follow the same procedure for the HEBE, we estimate , comparing this with the B-S relation, we find: sDF = 6(sDS-B) 2 the dimensionless sDF = 0.67, and DF = 1.33 Wm -2 K -1 , T4 = 2.23K (again not far from the data).We note that the ratio DF / DB-S ≈2 so that the estimates are close.
We can use this information to estimate lh in the HEBE.From the definition of DB-S as a thermal conduction coefficient per radian we obtain DB-S = K/R so that .To find the transport length, we can use , , to obtain: Alternatively, we can estimate lh from the global scale DH: We see that these lh estimates differ by a factor of bDB-S/DF ≈ b/2.Since typical numerical models with resolutions of hundreds of kilometers use kv ≈ 10 -4 m 2 /s, and kh ≈ 1m 2 /s, at least at these scales b ≈10 -2 so that the difference in the estimates may be large.For example since sDB-S ≈ 0.33, we find that the former yields lh ≈ 20km, while the latter yields, lh ≈ 4000 km.One way to reconcile the difference is to assume that b -that characterizes the horizontal-vertical effective diffusivity ratio -has a ( ) systematic scale dependence due to a difference in the scaling properties of kh and kv so that at global scales b ≈ 1 (this may arise as a consequence of the scaling anisotropic horizontal structure of the atmosphere at weather scales, notably of the horizontal wind field, the 23/9D model, [Schertzer and Lovejoy, 1985]).
A different (possibly additional) way of reconciling the estimates is to consider the potentially large (multifractal) intermittency of the diffusivities that introduces s strong scale effect.For example, [Havlin and Ben-Avraham, 1987], [Weissman, 1988], [Lovejoy et al., 1998]) show that in 1-D, the large scale effective thermal resistance rT -the inverse diffusivity -is the average of the small scale resistances.If we denote the spatial averages over a scale L by a subscript, and assume that the resistivity is scaling (scale invariant) up to planetary scales (denote this by R), then it will generally follow the following multifractal statistics: (40) Where the angle brackets denote statistical averages and Kr(q) is the moment scaling function that characterizes the scaling of the q th order statistical moment order of the thermal resistance.
The thermal resistance is proportional to the inverse thermal diffusivity, therefore the effective HEBE diffusive transport coefficient at scale L satisfies: (41) Finally, using we obtain: (42) Which relates the transport length at small scales L and planetary scales R. Depending on Kr(-1), the ratio can be quite small.For example, if the thermal resistivity statistics are taken as log-normal, then: so that so that .As discussed in appendix A, C1 ≈ 0.16 for the temperature in space (see also [Lovejoy, 2018]).Using this value as a guide, we find so that depending on the small scale resolution L, we can easily explain a factor of 10 or more increase in the effective transport length at large scales.Clearly the scale dependence of kh, kv is an important topic for future FEBE research.
3. The inhomogeneous heat equation

Babenko's method
The homogeneous heat equation in a semi-infinite domain is a classical problem and conductive -radiative surface boundary conditions naturally lead to fractional order operators, the HEBE and GHEBE.Although we have seen that fractional operators appear quite naturally, their advantages are much more compelling for the more realistic inhomogeneous equations relevant for the Earth.We therefore proceed to derive the inhomogeneous HEBE and GHEBE using Babenko's method.The more usual application is to find the surface heat flux given a solution to the conduction equation (see for example [Magin et al., 2004], [Chenkuan and Clarkson, 2018]), the following application appears to be original.
In the inhomogeneous case with t = t (x), lh = lh(x), lv = lv(x), a = a(x), there is no unique nondimensionalization.Therefore, we express the inhomogeneous anomaly heat equation with nondimensional operators as: (43) Where we have used and z is a time independent horizontal transport operator allowing for both advective and diffusive transport.Under the fairly general conditions, when z operates on the temperature field, it is proportional to the nondimensional divergence of the horizontal heat flux (discussed in part I, see eq. 4).Since the forcing is via the surface boundary condition rather than by an inhomogeneous term, eq.43 is mathematically homogeneous.
The first step in Babenko's method (see e.g.[Podlubny, 1999], [Magin et al., 2004]), is to factor the differential operator: (44) As usual, the general solution of a homogeneous equation is a linear combination of elementary solutions A+ and A-: The A+ solution leads to solutions that diverge at whereas A-leads to the required physical solutions with , ( [Podlubny, 1999]).Therefore we are interested in solutions to: putting z = 0 and using (part I, eq.22), we obtain: where Ts(t ,x) is the surface temperature anomaly and is the heat flux into the surface (the negative of which is the z component of the surface conductive (sensible) heat flux).Before interpreting the half order operator on the left, we can already give this equation a physical interpretation.When >0, sensible heat is forced into the Earth, some of it is stored in the subsurface (the term, the same horizontal position x but stored by heating up the subsurface, z<0), and some of the heat (the term), is transported horizontally to neighbouring regions (and conversely when <0).We can also understand the basic difference between the A+ and A-solutions: whereas the physically relevant A-solution correspond to energy storage and horizontal transport in the region z<0, the A+ solutions correspond to the region z>0 assumed to be devoid of conducting material.
The final step is to use the fact that the conductive heat flux is equal to the radiative imbalance (part I, fig.1): (48) Combining the equations 29, 30 we obtain the inhomogeneous Generalized Half-order Energy Balance Equation (GHEBE): (49) If needed, the internal field T(t,x;z), can be found by solving eq.49 for Ts(t, x) which is the z = 0 boundary condition for the full eq.43.We see that eq.49 reduces to the homogeneous GHEBE (eq.17) when t, lh, s, a are constant.
By comparing this derivation with that of the homogeneous GHEBE via the classical Laplace-Fourier transform method (section 2.1), it is clear that Babenko's method is very similar, but is more general.Whereas in the homogeneous equation, where the transforms reduce the derivative operations to algebra, the difficulty with Babenko's method is to find proper interpretations of the fractional operators.However, in the above, we assumed that t was only a function of position, so that ( ) = T t, x;0 ( ) Laplace (or Fourier) transform methods still apply in the time domain, in the next section we discuss the more challenging interpretation of the fractional inhomogeneous spatial operators.

The zeroth order high frequency GHEBE: the HEBE
Before discussing the inhomogeneous GHEBE, consider the case where the horizontal term lhz is small compared to ; below we argue that this is a good approximation for scales up to years and decades and greater than tens of kilometers (table 1, appendix A).Recall that the this horizontal transport term is in fact proportional to the divergence of the horizontal heat flux so that it may be small even when heat fluxes are significant [Trenberth et al., 2009].Alternatively, in globally averaged models, there are no horizontal inhomogeneities so that z = 0.In these cases ; and we obtain the inhomogeneous HEBE as a special case of the inhomogeneous FEBE: (50) We have written it with a general H since as in part I, an inhomogeneous version of the EBE may be obtained with H = 1.We have also used the Weyl derivative (i.e. from ) since this accommodated periodic or statistically stationary forcing as well as forcing starting at t = 0 (I this case we simply consider F = 0 for t≤0).Eq.50 shows that the HEBE only depends on the local climate sensitivity and the local relaxation time.We'll see below that explicit dependence on the horizontal transport (v, kh) and specific heat per volume rc is only important at scales somewhat smaller than the transport length scale (or alternatively at extremely long time scales, section 3.6).Before solving the HEBE, it is instructive to introduce the notation . is the equilibrium temperature that would be reached at time t, if at each location x, F was suddenly stopped and fixed at that value.With this notation, we may integrate both sides of eq.50 by order H, and multiply by t -H to obtain: (51) Written in this form, it is obvious that the temperature is constantly relaxing in a power law manner to (although if F is time dependent, equilibrium will in general never in fact be established).In the usual EBM special case (H = 1), the power law must be replaced by an exponential, the HEBE is obtained with H = 1/2.Since , the deviation from -the τ ∂ ∂t term (eq.50) -physically corresponds to the energy imbalance, as before, it is a power law, long memory energy storage term.
The FEBE is a linear differential equation that can be solved using Green's functions [Miller and Ross, 1993], [Podlubny, 1999].The solution is: where Gd,H is the H order Mittag-Leffler impulse response Green's function ( [Lovejoy, 2019a]).In general, Gd,H is only expressible in terms of infinite series, exceptions are the H = 1 EBE (Gd,1 = e -t ); and the H = ½ HEBE (eq.33).

Some features of stochastic forcing
The FEBE and the HEBE are examples of fractional relaxation equations; these have primarily been discussed in the context of deterministic forcings that start at t = 0.The corresponding stochastic fractional relaxation processes (in physics, "fractional Langevin equations", (FLE) see the references in [Lovejoy, 2019a]) -here corresponds to stochastic internal forcing.The FLE has received little attention, although [Kobelev and Romanov, 2000], [West et al., 2003] discuss the corresponding nonstationary random walks.The statistically stationary stochastic case that results when Weyl rather than Riemann-Liouville fractional derivatives are used is treated in [Lovejoy, 2019a], including the HEBE autocorrelation function and prediction problem (and its limits) when F is a Gaussian white noise.
To understand the noise driven HEBE, it is helpful to Fourier analyze it using [Lovejoy, 2019a], section 3.3 part I.At high frequencies, the derivative (energy storage) term dominates so that the temperature is a fractional integral (order H) of the forcing.At low frequencies, the derivative term can be neglected so that T ≈ sF implying that the equilibrium temperature follows the forcing and that s is indeed the usual climate sensitivity.
Alternatively, in real space, if F(t) is a unit step function Q(t) and s = 1, then for H ≠ 1 the long time relaxation to the equilibrium temperature response is a power law: (part I eq.33).Similarly, for small t, for H < 1, the impulse response is singular (part I eq.33).Due to this singularity, when F(t) is a Gaussian white noise, at high frequencies, T will be a fractional Gaussian noise (fGn) with exponent HfGn = H -½; averages over time Dt will behave as .When H ≤1/2 (HfGn≤0) and the resolution is increased ( ), this implies strong resolution dependencies (mathematically, small scale divergences) and so it is important in data analysis, including the estimation of the temperature of the Earth [Lovejoy, 2017].When forced by a white noise, the HEBE is exactly at the critical value HfGn = 0 corresponding to a "1/f" noise (research in progress indicates that it is at least close to a white noise).A particularly relevant aspect is that the correlation function and spectrum change very slowly from high to low frequencies [Lovejoy, 2019a].With data over a limited ranges of scales -e.g.months to decades -then, depending on the relaxation time t, the HEBE could mimic the FEBE with any H in the range 0< H ≤ ½ (hence -1/2≤ HfGn ≤0).It can therefore potentially account for the geographical variations in H reported in [Lovejoy et al., 2017] as being spurious consequences of geographical variations in t(x).
At global scales, the high and low frequency HEBE behaviours are close to observations.For example, the global value H = 0.5±0.2 was found for the long time behaviour needed to project the earth's temperature to 2100 [Hebert, 2017], [Hébert et al., 2020], and [Procyk et al., 2020] also using centennial scale global temperature estimates but using the FEBE directly, found the less uncertain H = 0.38±0.05;and using data at monthly and seasonal scales [Del Rio Amador and Lovejoy, 2019] found the value H = 0.42±0.03and used it).Appendix B discusses the spatial cross correlation matrix implied by the HEBE that is needed for example in calculating Empirical Orthogonal Functions (EOFs, or for the space-time macroweather model developped in [Del Rio Amador and Lovejoy, 2020b]).
Although the HEBE was derived for anomalies, these were not defined as small perturbations but rather as time-varying components of the full solution of the temperature (energy) equation with the time independent part corresponding to the climate state.The only point at which T was assumed to be small was with respect to the absolute local climate temperature about which the black body radiation was linearized, a fairly weak restriction on T. We could also mention that by allowing the albedo or other parameters to change in time, the HEBE could easily be extended to the study of past or future climates where it would broaden the spectrum potentially improving the modeling of glacial cycles.
An important feature of fractional differential operators is that they imply long memories, this is the source of the skill in macroweather forecasts ( [Lovejoy et al., 2015], [Del Rio Amador and Lovejoy, 2019]).The fractional term with the long memory corresponds to the energy storage process.In contrast, [Lionel et al., 2014] introduced a class of ad hoc Energy Balance Models with Memory (EBMM) whose (nonfractional) time derivative depends on integrals over the past state of the system.

The first order in space GHEBE
The HEBE is the GHEBE limit where horizontal transport effects are dominated by temporal relaxation processes and are ignored.Although this spatial scale depends on the time scale, appendix A estimates that at monthly time scales, this spatial scale is less ≈10 km and even at centennial scales it may only be only 100km or so.For these small spatial scales, we follow [Babenko, 1986], [Kulish and Lage, 2000], [Magin et al., 2004], and expand the square root operator using the binomial expansion: (for the expansion to be strictly valid, t must be a constant in time and in space; we have already assumed that is independent of time).As usual with Babenko's method, a rigorous mathematical justification is not available ( [Podlubny, 1999]), although recall that t, and lh are only functions of position so that for the temporal operator, Laplace and Fourier transforms techniques still work.
Considering the spatial part of the fractional operator, we see that it is weighted by the effective heat transport velocity V; as shown below, it plays the role of a small parameter (table 1, appendix A estimate it as ≈10 -4 m/s).Therefore, dropping the subscript "s" here and below, the GHEBE is: with the Weyl fractional derivatives (these are partial fractional derivatives).
Keeping only the spatial terms leading in the small parameter V, we have the first order (in space) GHEBE: (55) Or: (56) This equation is apparently similar to the usual transport equation.To see this, operate on both sides by , to obtain: Except for the factor ½, the half order derivative term and the "effective", (roughened) forcing, this is the usual transport equation.Nevertheless, although tempting, it would be wrong to think of this simply as a usual transport equation with an extra fractional term.The reason is that the extra term is not a small perturbation, it is dominant except at small spatial scales.
On the contrary, it is rather the classical transport terms that are small perturbations to the main HEBE.Alternatively, without the term, eq.59 is a generalized fractional diffusion equation (e.g.[Coffey et al., 2012]), although still with a key difference being that the fractional derivative is Weyl, not Riemann-Liouville (i.e. over the range to t, not 0 to t).

The equilibrium temperature distribution: the HEBE climate
The HEBE applies to time scales sufficiently short and to spatial scales sufficiently large that the horizontal temperature fluxes are too slow to be important, they are neglected.The first order correction (eqs.56, 57) makes a small improvement by giving a more realistic treatment of the small scale horizontal transport.However, a long time after performing a step increase of the forcing, the time derivatives vanish and a new climate state is reached.If the temperature followed the pure HEBE, the spatial equilibrium temperature distribution would be determined by setting the HEBE time derivative to zero: (58) Where the subscript "eq" indicates the long time equilibrium (climate) FEBE limit.However, appendix A shows thatdepending on the nature of the horizontal transport -at scales perhaps of the order of centuries, the horizontal heat fluxes will dominate the relaxation processes so that for very long times, this HEBE estimate is only approximate.

Equilibrium and approach to equilibrium in the inhomogeneous GHEBE
To understand the long time behaviour, we return to the GHEBE but perform a (long-time) binomial expansion of the halforder operator assuming that the transport terms dominate: (from here on we drop the "h" subscripts on l and the gradient operator).Again, to be strictly valid, t must be a constant so that and commute.We have to be careful since the advection length and relaxation times are functions of position (but not time) so that the spatial operators don't commute.Keeping terms to first order in time, we obtain: (60) To make progess, let's choose the transport operator so that its half powers are easy to interpret.The simplest approach is consider only diffusive transport and to use an isotropic fractional operator defined over the surface of the earth.For an arbitrary test function r, the corresponding order H fractional integral is: (for 0≤H≤d, where d is the dimension of space, here d = 2, see e.g.[Schertzer and Lovejoy, 1987], appendix A).This can be understood since in Fourier space, the Laplacian is and its inverse is , the "Poisson solver".
With the help of spherical harmonics, Appendix C generalizes the results of section 2.3 gives the corresponding operators and their fractional extensions on the surface of the sphere.
Applying eq.61 to the case d = 2 and H = 1 we have: Therefore, let us define a diffusive type transport operator and its inverse implicitly from its inverse half-order power: Hence let us define the half-order operator by: (64) With this definition the surface temperature equation 60 becomes: Where the range of the integration W = E is the entire surface of the earth.This equation has only superficial links to equations studied in the literature such as the "generalized fractional advection-dispersion equation" (e.g.[Meerschaert and Sikorskii, 2012], [Hilfer, 2000]).We can now consider the system reaching equilibrium after a step forcing F(x,t) = F0(x)Q(t), (increase by F0 (x) "turned on" at t = 0).At long enough times, the earth reaches equilibrium, the time derivative term vanishes and we obtain the equation for the equilibrium (climatological) temperatures: To obtain an approximate solution, let's now assume that differs from the climatological FEBE climate temperature Teq,FEBE(x) by a small perturbation dT(x). (67) then, using Teq(x) ≈ s (x)F0(x) in the integral, we obtain the approximation: is the slow, diffusive correction to the "instantaneous" (fast, high frequency), HEBE climate sensitivity s(x) that is estimated at usual (e.g.decadal) scales.As expected, since this is the long time solution after a step perturbation, it doesn't depend on t. lζ ( ) T eq x ( ) T eq x ( ) = T eq,HEBE x ( ) + δT x ( ) ; T eq,HEBE x ( ) = s x ( ) F 0 x ( ) Horizontal transport of heat redistributes the energy fluxes locally, but since the GHEBE is linear, it shouldn't affect the overall (global) energy balance.Let us check this by direct calculation of the globally averaged temperature.Averaging eq.66, we obtain: (69) Where the spatial averaging operator (overbar) is defined for an arbitrary function f.The average of the horizontal heat flux term yields: (70 Where KE is an unimportant constant from the x integration, independent of x'.The far right equality is an application of the divergence theorem on the surface E whose boundary is dE, ds is a vector parallel to the bounding line.But since the integration is over the whole earth surface (E), there is no boundary, hence the result.We conclude that while horizontal diffusion transports heat over the earth's surface, it does not affect the overall global radiation budget: .

Conclusions
Up until now, at macroweather and climate scales, the Earth's energy balance has been modelled using two classical approaches.On the one hand, Budyko -Sellers models assume the continuum mechanics heat equation, this yielding a 1-D latitudinally varying climate state.On the other hand, there are the zero-dimensional box models that combine Newton's law of cooling with the assumption of an instantaneous temperature-storage relationship.Both models avoid the critical conductive -radiative surface boundary condition; the former by ignoring heat storage, redirecting radiative imbalances meridionally away from the equator, the latter by postulating a surface heat flux that is not simultaneously consistent with the heat equation and energy conservation across a conducting and radiating surface (part I).
This two part paper re-examined the classical heat equation with classical semi-infinite geometry.In the horizontally homogeneous case (part I), the fundamental novelty is the treatment of the conductive -radiative boundary conditions, here (part II), it is the use of Babenko's method to extend this to the more realistic horizontally inhomogeneous problem.In both cases, the semi-infinite subsurface geometry is only important over a shallow layer of the order of the diffusion depth where most of the storage occurs (roughly estimated as ≈ 100m in the ocean, ≈<10m over land, see table 1 and appendix A).
T eq x ( ) ) T eq = T eq,HEBE The key result was obtained by using standard Laplace and Fourier techniques.It was shown quite generally that the surface temperatures and heat fluxes are related by a half-order derivative relationship.This means that if Budyko-Sellers models are right -that the continuum mechanics heat equation is a good approximation to the Earth averaged over a long enough timethat a consequence is that the energy stored is given by a power law convolution over its past history.This is a general consequence of the conductive -radiative surface boundary conditions in semi-infinite geometry and is very different from the box models that assume that the relationship between the temperature and heat storage is instantaneous.Although the system itself is classical, this result may be viewed as a nonclassical example of the Mori-Zwanzig mechanism in which system parameters that are not modelled explicitly (here, the subsurface temperatures) imply long (power law) memories for the modelled parameters (here, the surface temperatures).This is in contrast to conventional short (exponential) memory assumption.It implies that any part of the Earth system that exchanges energy both radiatively and conductively into a surface should be modelled with fractional rather than integer ordered derivatives.A far reaching consequence is that classical dynamical systems approaches based on integer ordered differential equations are not necessarily pertinent to the climate system.
If we ignore horizontal heat transport (part I), an immediate consequence of half order storage is that the temperature obeys the Half-order Energy Balance Equation (HEBE) rather than the classical first order EBE.Depending on the space-time statistics of the anomaly forcing, the HEBE justifies the current Fractional EBE (FEBE) based macroweather (monthly, seasonal) temperature forecasts [Lovejoy et al., 2015], [Del Rio Amador and Lovejoy, 2019], [Del Rio Amador and Lovejoy, 2020a;Del Rio Amador and Lovejoy, 2020b] that are effectively high frequency approximations to the FEBE).Similarly, the low frequency (asymptotic) power law part can produce climate projections with significantly lower uncertainties than current GCM based alternatives ( [Hebert, 2017], [Hébert et al., 2020] and work in progress directly using the HEBE, [Procyk et al., 2020]).
When the system is periodically forced, the response is shifted in phase -and borrowing from the engineering literature -the surface is characterized by a complex thermal impedance that we showed is equal to the (complex) climate sensitivity.In part I, we gave evidence that this quantitatively explains the phase lag (typically of about 25 days) between the annual solar forcing and temperature response.
In this second part, we investigated the consequences of horizontal heat transport, first in a homogeneous medium with inhomogeneous forcing first on a plane and then -permitting a direct comparison with the usual Budyko-Sellers approachon the sphere (section 2).In section 3 we considered inhomogeneous material properties (including variable diffusion lengths, relaxation times, and climate sensitivities).While Laplace and Fourier techniques can still be used in time, they cannot be used in space.However, the extension to inhomogeneous media was nevertheless possible thanks to Babenko's powerful (but less rigorous) operator method.Whereas in part I, the homogeneous fractional space-time operator was given a precise meaning, here -following Babenko -the corresponding inhomogeneous operator was interpreted using binomial expansions for both the short and long time limits and yield 2D energy balance models.Part II thus allows us for the first time to extend energy balance models to 2-D, allowing the treatment of regional temporal anomalies.
The expansions depend both on the space and time scale and on a dimensional parameter: the typical horizontal transport speed (V), estimated as ≈ 10 -4 m/s (appendix A).The zeroth order expansion in time limit yielded the inhomogeneous HEBE, the first order correction yielded an equation that superficially resembled the usual heat equation but instead had a leading halforder time derivative term.Based on the analysis of NCEP reanalyses (appendix A), it was argued that at spatial scales larger than hundreds of kilometers, that these approximations are likely to be useful for years, decades, and perhaps longer.However, for studying climate states -defined for example as the equilibrium state for forcings that are increased everywhere in step function fashion -we required low, not high frequency expansions and these are based on fractional spatial operators.We defined inhomogeneous fractional diffusion operators in both flat space and on the sphere (appendix C), and derived equations for both the equilibrium limit and the approach to the limit.We showed that (as expected) they conserved energy and that the low frequency climate sensitivity is somewhat different from that estimated at higher frequencies (from the EBE or HEBE).
The EBE and HEBE are the H = 1, H = 1/2 special cases of the Fractional EBE (FEBE) that was recently introduced as a phenomenological model [Lovejoy et al., 2020] (see also [Lovejoy, 2019a], [Lovejoy, 2019b]) with empirical estimates H ≈ 0.4 -0.5, i.e. very close to the HEBE.Although only a special case, the HEBE illustrates the general features of the FEBE fractional-order energy storage term and power law long memories, in [Lovejoy, 2019a] discussed the statistical properties of the FEBE driven by Gaussian white noise (a model for the internal variability forcing) showing that the high frequency limit is a process called fractional Gaussian noise (fGn).In the special HEBE case with H = 1/2, the fGn temperature response has exactly a high frequency 1/f spectrum that is cut-off at the relaxation time (empirically of the order of a few years).[Lovejoy, 2019a] developed optimal predictors and determined the predictability skill.
Whereas the more general FEBE is essentially a phenomenological model up until now justified by the hypothesized scale invariance of the energy storage mechanisms ( [Lovejoy et al., 2020]), the HEBE follows directly and quite generally from the continuum mechanics heat equation, thus giving it a more solid theoretical basis.However, the work here suggests another way to obtain the FEBE: to replace the classical heat equation by its fractional generalization, the fractional heat equation, a possibility that we explore elsewhere.
As a final comment, we should mention that although this paper focused on the time varying anomalies with respect to a time independent climate state, our approach opens the door to new methods for determining full 2-D climate states (generalizations of the 1-D Budyko-Sellers type climates) but also to determining past and future climates and the transitions between them.This is because the definition of temperature "anomalies" is very flexible.For example, we could first apply the method to determining the existing climate by fixing the forcing at current values and solving the time independent transport equations.Then, the long term effect of changes such as step function increases in forcing could be determined from the GHEBE anomaly equation (section 3.5) which regionally corrects the local climate sensitivities for (slow) horizontal energy transport effects.Nonlinear effects that can be modelled by temperature dependent forcings (i.e. ) can easily be introduced.Other nonlinear effects needed to account for Milankovitch cycles could thus easily be made, the primary difference being the half-order derivatives and the scaling that they imply.Indeed, the In the macroweather regime, the temporal temperature fluctuation at time scale Dt is where is the anomaly averaged over scale Dt; empirically this is valid over the macroweather regime i.e. up to 10 -30 years in the industrial epoch (see e.g.[Lovejoy and Schertzer, 2013], [Lovejoy, 2013], [Lovejoy et al., 2017]).The typical fluctuation can be estimated by the RMS anomaly: Where the overbar is the average over all the anomalies in a time series at a single location x.Dt1 is a convenient reference time, here taken as 1 month.Empirically, the exponent Ht ≈ 0 to -0.2; this is similar to the high frequency result Ht = 0 (i.e. for Dt<t) predicted from the HEBE with white noise forcing, valid for Dt <t.Hence for our present purposes the typical time derivative is: (74) This is the resolution Dt time derivative.Since typical north-south gradients are larger than typical east-west ones, the meridional (y) component of the transport is dominant, so that we will focus on it: (75) Hence the meirdional contributions to the ratios ra, rd are: Where , is the relative fluctuation in the RMS temperature at time scale Dt, spatial scale Dy and -since we are only interested in an order of magnitude -we took a ≈ ay.The estimate of the diffusive term uses a finite difference approximation to the Laplacian.lh is horizontal diffusion length and a is the nondimensional advection speed v/V (V = lh/t, see below).To gauge the order of magnitudes, in the far right term of eq.76, we took the absolute value so that the 625 result is an upper bound.
Climate sensitivity s water ≈ 4x10 Table 1 summarizes the dimensional and nondimensional parameter estimates, the final step is to estimate values of the gradient and Laplacian terms (eq.76).Since s -and hence log s -are the amplitudes of temporal noises; these amplitudes vary 630 stochastically from one spatial location to another.Due to the space-time scaling of the temperature anomalies (analysed in [Lovejoy and Schertzer, 2013]), we expect that the statistics of the logarithms (eq.76) to follow power laws up to large scales.To quantify this, we used NCEP reanalysis data at 2.5 o resolution from 1948 to present, and after removing the low frequency anthropogenic trend, we estimated the RMS temperature anomalies at each pixel; s(x).In fig.6, we then calculated spatial zonal and meridional fluctuations Dlogs(Dx), Dlogs(Dy), and from these their root mean square (RMS) values.From the figure, we see that to a good approximation: (77) The fluctuations we used are Haar fluctuations, but because Hx ≈ Hy > 0, they are nearly equal to difference fluctuations [Lovejoy and Schertzer, 2012].We see that the zonal and meridional lines are roughly parallel: with a "trivial" horizontal anisotropy factor ≈ 5 (typical north-south fluctuations are 5 times larger than typical east-west ones).Although, H = 1/2 is the value corresponding to Brownian motion, the actual variability is highly intermittent (spiky), so that unlike the temporal fluctuations, these spatial increments are far from Gaussian; it is not Brownian motion.Multifractal analysis indicates that the intermittency parameter (the codimension of the mean) C1 ≈ 0.16 which is very high, reflecting the strong spatial fluctuations as we move from one climate zone to another [Lovejoy and Schertzer, 2013], [Lovejoy, 2018], [Lovejoy, 2019b].
Since the north-south gradients are much stronger than the east-west ones, we can estimate the gradients and Laplacians by using the y direction fluctuations: at scale Dy: (78) (79) Since LNS ≈ 3x10 6 m , over most of the range of Dy, so that the ratio of advection to diffusion is so that advection dominates diffusion for .Taking a ≈ 1, it is dominant for .

Fig. 4 :Fig. 6 :
Fig. 4: This is the step response in time and (circular) step in space for conductive-radiative forcing.Lines for t = 0.01 (bottom), 0.2, 0.4, ... 1.6 (black, bottom to top, the thick black line is for ( equilibrium).The nondimensional forcing is the rectangle (from unit circular forcing).Also shown (top dashed) is the equilibrium when the forcing is purely due to unit conductive heating over the unit circle.

Table 1 :
Parameter estimates from part 1 section 3.1.2,see section 2.3 for some planetary scale estimates.