{'fraction_non_alphanumeric': 0.03982136868303532, 'fraction_numerical': 0.017721921169223165, 'mean_word_length': 4.783380971825235, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The half-order energy balance equation – Part 2: The inhomogeneous HEBE and 2D energy balance models

In Part 1, I considered the zero-dimensional heat equation, showing quite generally that conductive– radiative surface boundary conditions lead to half-ordered derivative relationships between surface heat fluxes and temperatures: the half-ordered energy balance equation (HEBE). The real Earth, even when averaged in time over the weather scales (up to ≈ 10 d), is highly heterogeneous. In this Part 2, the treatment is extended to the horizontal direction. I first consider a homogeneous Earth but with spatially varying forcing on both a plane and on the sphere: the new equations are compared with the canonical 1D Budyko–Sellers equations. Using Laplace and Fourier techniques, I derive the generalized HEBE (the GHEBE) based on half-ordered space–time operators. I analytically solve the homogeneous GHEBE and show how these operators can be given precise interpretations. I then consider the full inhomogeneous problem with horizontally varying diffusivities, thermal capacities, climate sensitivities, and forcings. For this I use Babenko’s operator method, which generalizes Laplace and Fourier methods. By expanding the inhomogeneous space–time operator at both high and low frequencies, I derive 2D energy balance equations that can be used for macroweather forecasting, climate projections, and studying the approach to new (equilibrium) climate states when the forcings are all increased and held constant.

The homogeneous case is quite classical and can be treated with standard Laplace and Fourier techniques, it leads to the 40 (horizontally) Generalized HEBE: the GHEBE. Although the GHEBE has a more complex fractional derivative operator, 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. This done by using of Babenko's method [Babenko, 1986] which is essentially a generalization of the Laplace and Fourier transform techniques. The challenge with Babenko's method is to interpret the inhomogenous fractional operators. Following Babenko, we do this using both high and low 45 frequency expansions corresponding respectively to processes dominated by storage and by horizontal heat transport. The long time limit describes the new climate state that results when the forcing is increased everywhere and held fixed. We also include several appendices focused on empirical parameter estimates (appendix A), the implications for two point and spacetime temperature statistics (when the system is stochastically forced, internal variability, appendices B, C), and finally (appendix D), the changes needed to account for the Earth's spherical geometry, including the definition of fractional operators 50 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): 55 (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: 60 (2) r, c are the fluid densities and specific heats l is the climate sensitivity and F is the anomaly forcing. The initial conditions are T = 0 at (all t), and (Riemann-Liouville) or below, . In part I, we nondimensionalized the zero-dimensional homogeneous operators by nondimensionalizing time by the relaxation 65 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 l = 1 (equivalent to using dimensions of temperature for the forcing F), we obtain: 70 (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) 75 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 D, we consider the changes needed to account for spherical geometry.
When , the solution and where Gd is the impulse 80 (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 85 (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: 90 (7) 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: 95 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. 30, part I. The 100 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 center of the impulse response with nondimensional velocity a. As 105 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. 31, part I). From these, the general surface results including advection are obtained with , i.e. 110 . 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 115 power law approach to thermodynamic 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" thermodynamic equilibrium response 120 discussed below (eq. 20). Note that the leading term in is independent of r, and the leading term in the approach to thermodynamic 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 125 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 130 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: ("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 140 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.

Thermodynamic equilibrium
If F(t) = 0, then the system is at equilibrium and will stay there. However, if F is a step function in time, then as , a 145 new equilibrium will be established. At equilibrium, d/dt = 0, so that the conjugate variable p = 0. With this and a = 0 in eq. 15, we obtain the equation for the (spatial) surface impulse response for thermodynamic equilibrium (subscript "therm"): i.e. the same as eq. 4 but with p = 0 (and a = 0) hence: 150 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: 155 ; 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) 160 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 thermodynamic 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: 165 (22) P1(r) is the corresponding indicator function. We now use the transform pair to perform the convolution: (J1 is the first order Bessel function of first kind). Taking the limit we obtain the thermodynamic equilibrium temperature distribution. Alternatively we could find it directly by from eq. 19: 170 (24) 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 thermodynamic equilibrium when the forcing is purely due to unit conductive heating over the unit circle. The difference between the dashed and the thick thermodynamic equilibrium curves are purely 175 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 thermodynamic equilibrium. This explains why Brunt instead used temperature forcing boundary conditions. Here, in two horizontal dimensions, boundary conditions that impose a fixed G therm,δ r;0 temperature over the circle are problematic since they imply infinite horizontal temperature gradients and infinite horizontal heat fluxes). 180 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).

Babenko's method 185
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 naturally, their advantages are much more compelling for the more realistic inhomogeneous equations relevant for the Earth.
We will therefore now 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], 190 [Chenkuan andClarkson, 2018]), the following 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: Where we have used and z is a time independent horizontal transport operator allowing for 195 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. 25 is mathematically homogeneous.

(27)
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: (28) 205 putting z = 0 and using eq. 26, 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 210 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. 215 The final step is to use the fact that the conductive heat flux is equal to the radiative imbalance (part I, fig. 1): Combining the equations 29, 30 we obtain the inhomogeneous Generalized Half-order Energy Balance Equation (GHEBE): If needed, the internal field T(t,x;z), can be found by solving eq. 31 for Ts(t, x) which is the z = 0 boundary condition for the 220 full eq. 25. We see that eq. 31 reduces to the homogeneous GHEBE (eq. 17) when t, lh, l, 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 225 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  230 1, appendix A). 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: 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 235 stationary forcing as well as forcing starting at t = 0 (I this case we simply consider F = 0 for t≤0). Eq. 32 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.5). Before solving the HEBE, it is instructive to introduce the notation . is the equilibrium temperature that 240 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. 32 by order H, and multiply by t -H to obtain: Written in this form, it is obvious that the temperature is constantly relaxing in a power law manner to (although if F and is time dependent, equilibrium will in general never in fact be established). In the usual EBM special case (H = 1), the power 245 law must be replaced by an exponential, the HEBE is obtained with H = 1/2. Since , physically the deviation from -the term (eq. 32) -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: 250 where G0,H is the H order Mittag-Leffler impulse response Green's function ( [Lovejoy, 2019a]). In general, G0,H is only expressible in terms of infinite series, exceptions are the EBE (G0,1 = e -t ); and the HEBE with in the notation above (eq. 31, part I).
The corresponding step response G1,1/2 = GQ is the integral of G0,1/2 (eq. 32, part I), it describes relaxation to thermodynamic 255 equilibrium when F is a step function; similarly, the ramp (linear forcing) response G2,1/2 (eq. 32, part I), is the integral of the step response.

Some features of stochastic forcing
The FEBE and the HEBE are examples of fractional relaxation equations; these have primarily been discussed in the context 260 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 have 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 265 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 and appendix C. 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 ≈ lF implying that the equilibrium temperature follows the forcing and that l is indeed the usual climate sensitivity. 270 Alternatively, in real space, if F(t) is a unit step function Q(t) and l = 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) this implies strong resolution dependencies (mathematically, 275 divergences) when the resolution is increased ( ) and so 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 (note that the Earth's internal variability forcing is not necessarily a white noise, it might have a different scaling behaviour). 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, 280 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  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]; the value H 285 = 0.42±0.03 for the internal macroweather variability needed to make monthly and seasonal forecasts [Del Rio Amador and  (note that this was inferred by make the usual assumption that the internal forcing F is a Gaussian white noise, and this may not be the case). Appendix B discusses the spatial cross correlation matrix implied by the HEBE that is needed for example in calculating Empirical Orthogonal Functions (EOFs).
We could also mention that if F is spatially statistically homogeneous and independent of the parameters l, t, then not only 290 will the macroweather temperature fluctuations be well reproduced, but also, up to the relaxation time, the temperature may easily respect a space-time symmetry called space-time statistical factorization, ("STSF"; e.g. where R represents the autocorrelation function), see appendix C. Empirically, the STSF is at least approximately obeyed by space-time temperature and precipitation fluctuations ( [Lovejoy and de Lima, 2015]), and if respected, the STSF has important implications for macroweather temperature forecasting. 295 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 300 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 , [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 305 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 310 scale is of the order of ≈1 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: 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 320 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: (37) 325 Or: This equation is apparently similar to the usual transport equation. To see this, operate on both sides by , to obtain: (39) 330 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 335 the term, eq. 41 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 HEBE thermodynamic equilibrium
The HEBE applies to time scales sufficiently short and to spatial scales sufficiently large that the horizontal temperature fluxes 340 are too slow to be important, they are neglected. The first order correction (eqs. 38, 39) 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 pattern for thermodynamic equilibrium would be determined by setting the HEBE time derivative to zero: Where the subscript "c" indicates the long time (climate) FEBE limit. However, appendix A shows that -depending on the nature of the horizontal transport -at scales perhaps of the order of millennia, 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 350
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). 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 355 first order in time, we obtain: 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: 360 (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".
Note that eqs. 42, 43 involve ½ order inverse Laplacians which are H = 1 (rather than H = ½) isotropic integrals (eq. 43). Applying eq. 43 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: 370 (45) Hence let us define the half-order operator by: With this definition the surface temperature equation 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 thermodynamic equilibrium and the time derivative term vanishes and we obtain the equation for the climatological temperatures: 380 To obtain an approximate solution, let's now assume that differs from the climatological FEBE climate temperature Tc,FEBE(x) by a small perturbation dT(x).
in the integral, we obtain the approximation: 385 (50) is the slow, diffusive correction to the "instantaneous" (fast, high frequency), HEBE climate sensitivity l(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.
Horizontal transport of heat redistributes the energy fluxes locally, but since the GHEBE is linear, it shouldn't affect the overall 390 (global) energy balance. Let us check this by direct calculation of the globally averaged temperature. Averaging eq. 48, we obtain: Where the spatial averaging operator (overbar) is defined for an arbitrary function f. The average of the horizontal heat flux term yields: 395 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: . 400

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 holds, this yields 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 405 conductive -radiative surface boundary conditions; 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 novelty is the treatment of the conductive -radiative boundary conditions, here (part II), it is 410 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).
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 415 right -that the continuum mechanics heat equation is a good approximation to the Earth averaged over a long enough timethen 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 condition and is very different from the box models that assume that the relationship between the temperature and heat storage is instantaneous.
If we ignore horizontal heat transport (part I), an immediate consequence of half order storage is that the temperature obeys 420 the Half-order Energy Balance Equation (HEBE) rather than the classical order one EBE. The implied long time storage behaviour explains the success of scaling based climate projections and, the implied short time behaviour potentially explains the success of macroweather forecasts that exploit it. 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 425 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 (section 2) and then, more generally with inhomogeneous material properties (including variable diffusion lengths, relaxation times, and climate sensitivities, section 3). Laplace and Fourier techniques and not so useful here, but the extension to inhomogeneous media was nevertheless possible thanks to Babenko's powerful (but less rigorous) operator 430 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.
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 435 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 thermodynamic 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 440 operators. We defined inhomogeneous fractional diffusion operators in both flat space and on the sphere (appendix D), and derived equations for both the thermodynamic 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 445 phenomenological model [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. [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 450 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 unitl now justified by the hypothesized scale invariance of the energy storage mechanisms, 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: 455 to replace the classical heat equation by its fractional generalization, a possibility that we explore elsewhere. Part II allowed us for the first time to extend energy balance models to 2-D, allowing the treatment of regional temporal anomalies. Depending on the space-time statistics of the anomaly forcing, the HEBE justifies the current Fractional EBE (FEBE) based macroweather (monthly, seasonal) temperature forecasts , [Del Rio Amador and Lovejoy, 2019]. Similarly, the low frequency (asymptotic) power law part can produce climate projections with significantly lower uncertainties than current 460 GCM based alternatives ( [Hebert, 2017] and work in progress directly using the HEBE with R. Procyk).
This work was performed in the spirit of Budyko-Sellers models in which the Earth system is averaged over scales longer than typical lifetimes of planetary scale weather structures. Following Budyko-Sellers, the key physical assumption was that the resulting averaged system is a continuum system, thus justifying use of the general continuum mechanics heat equation. From this, the GHEBE and HEBE follow from the surface conductive-radiative boundary condition. In as much as GCMs (that are 465 based on continuum mechanics) reproduce the same statistics as the noise -or anthropogenically forced FEBE and HEBE, the continuum hypothesis is plausible.
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 470 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. 475 ) 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 power law relaxation processes implied by the GHEBE suggests straightforward explanations for the observed power law climate regime spanning the range from centennial to Milankovitch scales.

Acknowledgements 480
I acknowledge discussions with L. Del Rio Amador, R. Procyk, R. Hébert, D. Clarke and C. Penland. This is a contribution to fundamental science; it was unfunded and there were no conflicts of interest. The data used in appendix A are from the NOAA website: https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html.

Appendix A: Empirical analysis of the horizontal structure
In order to apply our results to the Earth, we need some idea of the magnitudes of various terms in our equations. To 485 start with, recall that our model is of the Earth system at macroweather and climate time scales i.e. all relevant quantities are averaged over the weather scales ≈ 10 days or longer. The resulting averaged system is then treated as a continuum and the general continuum mechanics heat equation is applied. In this, we essentially follow the Budyko-Sellers approach and consider that the diffusive transport is characterized by eddy (not molecular) diffusivities and that the vertical structure of this averaged continuum is homogeneous (although it may vary considerably from place to place in the horizontal). Unlike Budyko-Sellers 490 that treat the vertical as negligibly thick -they don't consider it at all -our key difference is that we assume that it has a thickness of the order of a few diffusion depths, and then we apply the key conductive-radiative surface boundary condition.
Probably the most important aspect is to estimate the relative importance of the temporal relaxation (and storage) terms in comparison to the horizontal transport terms (see eq. 25). Indeed, for judging their relative importance, the key parameter is the ratio of the transport to relaxation terms r: 495 Where a is the magnitude of the dimensionless advection velocity vector a = v/V. When r≪1, the transport term is small compared to the temporal term, conversely when r≫1. In order to quantify this, it is convenient to consider the advective ("a") and diffusive ("d") terms as well as their derivatives individually: (54) 500 In the macroweather regime, the temporal temperature fluctuation at time scale Dt is where is the anomaly average 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, 2013], ). The typical fluctuation can be estimated by the RMS anomaly: (55) 505 Where the overbar is the average over all the anomalies in a time series at a single location x. Dt1 is convenient reference time, here taken as 1 month. Empirically, the exponent Ht ≈ 0 to -0.2; this 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: This is the resolution Dt time derivative. Since typical north-south gradients are larger than typical east-west ones, the meridional (y) component is dominant, so that we will focus on it: Hence: (58) 515 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 anomaly relaxation 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 of eq. 58, we suppressed the signs.
Alternatively we can use kv = t/(rcl) 2 and the global estimates of t ≈ 10 8 s ( [Hebert, 2017], work in progress with R. Procyk, or part I, section 3.3). From these, we obtain kv ≈ 10 -5 m 2 /s which is close to the model values. In conclusion, using kv ≈ 10 -5 -10 -4 m 2 /s yields lv ≈ 30 -100m, lh ≈ 10 km. Consequently, the diffusive based velocity parameter is V ≈ lh/t ≈ 10 -4 m/s. 530 The best transport model -diffusive, advective -or both -is not clear, therefore let us estimate the magnitude of the advective velocity v assuming that it dominates the transport. The appropriate value is not obvious since most models just use eddy diffusivity -not advection -for transport. One way -for example [Warren and Schneider, 1979] -is to note that typical meridional heat fluxes are of the order of 100 W/m 2 over meridional bands whose temperature gradients DT are several degrees If this heat is transported by advection, it implies v ≈ Qa/(rcDT) ≈ 10 -5 -10 -4 m/s (eq. 4, part I), hence, using V ≈ 10 -4 m/s 535 (above), we find a = v/V ≈ 0.1 -1.
With these dimensional and nondimensional parameter estimates, the final step is to estimate values of the gradient and Laplacian terms (eq. 58). Since s -and hence log s -are the amplitudes of temporal noises; these amplitudes vary stochastically from one spatial location to another. Due to the spatial scaling, we expect that their statistics 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 540 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: The fluctuations are Haar fluctuations, but because Hx ≈ Hy > 0, they are nearly equal to difference fluctuations [Lovejoy and 545 Schertzer, 2012]. We see that the zonal and meridional lines are roughly parallel: with a "trivial" horizontal anisotropy factor ≈ 5. Although, H = 1/2 is the value corresponding to Brownian motion, the actual variability is highly intermittent (spiky), so that the 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, 2018], [Lovejoy, 2019b]. 550 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: Since LNS ≈ 3x10 6 m , over most of the range of Dy, so that the ratio of advection to diffusion is 555 so that advection dominates diffusion for . Taking a ≈ 1, it is dominant for .
Using lh ≈10 4 m, LNS ≈ 3x10 6 m, Hy = 1/2, V =10 -4 m/s we find approximately critical length scales that yields unit ratios: Where Dt is measured in seconds, Dy in meters. When the typical distances exceed these critical distances (i.e. when Dy>Dyc), 560 we have r<1 so that the temporal derivative terms dominate over the horizontal transport. For Dt = 1 month, we have Dyc,a ≈ 0.1m, and Dyc,d ≈ 200m, so that unless the distances are very small, the temporal (storage) terms are indeed dominant. Even over much longer time scales -e.g. Dt ≈30 years (10 10 s), they dominate for distances greater than ≈Dyc,a ≈ Dyc,d ≈ 10 km.
Alternatively, we could estimate the time scales needed so that the critical transport scale is 1000km. From the same equations, we obtain estimates of 300 years (advection), 30,000 years (diffusion). Note however that in the anthropocene, for periods Dt 565 ≈> 10 years, that the temporal fluctuations start to grow (i.e. the empirical relations eqs. 60, 61 will break down); nevertheless, the above scaling relations for the internal variability may hold to much longer times .
In summary, from eq. 62, we conclude that for the larger scales >>≈10 km, that r≪1 and that the HEBE may apply except for time scales ≫t: the only explicit role of kh, kv, r, c is to determine the limits of validity of the HEBE via lh, a. When the HEBE is valid, only the relaxation time t and the climate sensitivity l are relevant. 570

Appendix B: The HEBE cross-correlations
The temperature anomaly cross-correlation function (a matrix when the temperature is discretized on a grid), is commonly used in climate science, notably to determine Empirical Orthogonal Functions (EOFs). These can be determined from the HEBE (or GHEBE if needed) once a forcing model is given. Let us first consider that the climate sensitivities and relaxation times are deterministic characterizations of the local properties at points x1, x2. In this case, for the HEBE, any correlations 575 between the temperature anomalies at those points will arise because of correlations in the forcing F(x,t). We now consider simple deterministic and stochastic forcings.

a) Deterministic forcing:
The simplest model is to take complete spatial correlation correlations, with a step function (Q(t)) forcing at t = 0, but different at each position x: 580 The temporally averaged cross-correlation can be determined by: Recalling that G1,1/2 (= GQ) is the step response, the integral of G0,1/2 (= Gd) since we have: Hence: (66)

b) Stochastic forcing: 590
A convenient model of pure internal variability, is to assume that the forcing is statistically stationary in time with the following forcing cross-correlations: (67) (the "<.>" symbol indicates ensemble, statistical averaging). This implies a stationary temperature cross-correlation: Note the general symmetry property ; we only need to determine R for Dt>0. For statistically stationary forcing, is the anomaly cross-correlation needed -for example -for constructing Empirical Orthogonal Functions (EOFs).
The easiest way to relate RF and RT is via their spectra. Let us define the transform pairs: (69) 600 similarly for the forcing F (the circonflex indicates Fourier Transform). Then: (70) (this is true for the Weyl fractional derivatives used here, [Podlubny, 1999]). So that the impulse response is: The solution to the HEBE at two different points x1, x2 is: Where the asterix indicates complex conjugate. Multiplying and taking ensemble averages and assuming that the forcingand hence responses -are statistical stationary, we obtain: 610 (73) Where: Therefore: (75) 615 A special case that is useful later, is when x1 = x2 = x, which yields the spectrum ET at the point x: Using a partial fraction expansion of eq. 75, we obtain: 620 By inverting the Fourier transform, this can be used to determine the real space transfer function . Using contour integration, it is convenient to convert the inverse Fourier transforms into Laplace transforms for Dt > 0:

(78)
For Dt<0, use . The spatial cross-correlation, temporal autocorrelation function of the temperature is therefore: Where the "*" indicates convolution. 630 The basic Laplace transforms in eq. 78 can be expressed in terms of higher mathematical functions as follows (all for t>0): (80)

635
The ip comes from integrating half way around the pole at the origin. Note that both the Exponential Integral (EI) and the incomplete Gamma functions have log divergences at the origin. If needed, these formulae can be combined to obtain a complete analytic expression for , which can then be used to determine the temperature correlations if the forcing correlations are known: where the asterix is the temporal convolution. 640 The special case x1 = x2 i.e. with t1 = t2 = t, is a little simpler: Whose Fourier transform is: Evaluating the integral for g(Dt) using the Laplace transform formulae (eq. 80): 645 Note the small scale log divergence, this is important when the forcing is a white noise, see [Lovejoy, 2019a]. The temporal 650 autocorrelation at the point x is thus: However, in general, the Fourier relations are easier to deal with.

Appendix C: Statistical Space-Time Factorization
At high frequencies (i.e. Dt<t), and empirically over the macroweather regime up to a decade or more ( [Lovejoy and de Lima,655 2015]), both precipitation and temperature anomalies (at least approximately) respect a space-time symmetry called "spacetime statistical factorization" ("STSF").
For example, for the autocorrelation function R, this implies . If obeyed, this factorization implies important simplications in regional macroweather forecasting: it is therefore interesting to investigate the implications HEBE for the STSF hypothesis.
The easiest way to approach the STSF is to consider that the forcing and relaxation times t(x) and sensitivities l(x) are 660 stochastic fields that are statistically homogeneous in space so that the correlation functions can be written: . If we assume that the forcing is statistically independent of the temperature, then, taking the high frequency limit of in eq. 75: we obtain: 665 From this, see that if the forcing factorizes then the temperature autocorrelation function also factorizes: Where is the autocorrelation function of lt -1/2 (the term in square brackets in eq. 87). From here, the inverse Fourier transform of and gives the real space version of the STSF symmetry. Notice that at the STSF hinges on the factorization approximation for and at low w, it breaks down.

Appendix D: Fractional Integration on the sphere
At long enough time scales, the spatial transport of heat is important and the spherical geometry of the Earth must be taken 675 into account. The standard way (see e.g. the review [North et al., 1981]) is to use spherical harmonics. In Appendix 5D of  these were used to define fractional integrals on the sphere, necessary in order to produce the corresponding multifractal cloud and topography models (see also [Landais et al., 2019]). Spherical harmonics are particularly convenient when the heat transport is diffusive, involving fractional Laplacians. In section 3.5.2, these were defined in real space by taking the domain of integration to be a sphere. In this appendix we discuss an alternative method of 680 spherical fractional integration that may have theoretical and practical advantages.
The Laplacian on a sphere ( ) is the angular part of the Laplacian in spherical coordinates, it is obtained by expressing the Laplacian in spherical coordinates and setting the radial derivatives to zero: The spherical harmonics form a complete orthogonal basis, so that any function on the sphere can be uniquely expressed in terms of a spherical harmonic expansion: Where the are the coefficients of the expansion without fractional integration (i.e. of order 0, indicated in the superscript). 695 This suggests the following definition for a fractional spherical integration order H of a spherical harmonic: , for the HEBE, we take H = 1 which corresponds to the ½ power of the inverse Laplacian. We have excluded the value n = 0 since when H>0, the filter divergences; since , this component corresponds to the mean. Therefore the above definition is adequate for mean zero anomalies. Alternatively, the mean can be removed and 700 taken care of separately, see below. With this definition, the fractional integral of the zero mean function f is: i.e. a filter in spherical harmonic space, analogous to the Fourier filter for an isotropic fractional integration in Cartesian coordinates.
The definition of the fractional Laplacian (eq. 93, 94) is adequate when the horizontal transport coefficients are constant, 705 but in section 3.5, we saw that more generally, the half order divergence operator was written: i.e. there was an extra multiplication by the spatially varying diffusion length . In flat (Cartesian) coordinates, such real space multiplications correspond to Fourier space convolutions so that this operator can also be conveniently expressed in Fourier space. However, with spherical harmonics, this simplicity is lost: although isotropic real space convolutions can still be performed by filtering the harmonics, real space multiplications no longer correspond to convolutions of harmonic coefficients, 710 the closest spherical harmonic equivalent is much more complicated, it involves Clebsch-Gordon coefficients.
A method of fractionally integrating the mean (n = 0) component was developed for the purpose of multifractal modeling in Appendix 5D of . There, a different definition of fractional integrals on the sphere was proposed: a convolution with the function Q -(2-H) , where Q is the angle between two points subtended at the center of the sphere. The function Q -(2-H) / G(H/2) was numerically expanded in spherical harmonics and the convolution was again 715 performed by filtering the coefficients (the constant G(H/2) is needed so that the normalization is the same as for the definition eq. 92). The main difference between the two definitions is that the latter can be directly applied to fields with nonzero means.
With definition, the H order fractional integral of a constant function on the sphere (representing the nonzero mean), is simply the value multiplied by which for the HEBE H =1 case, reduces to (1/2) 1/2 Si(2p) where Si is the standard sine integral function. However for the coefficients n≥1, numerical tests show that the two definitions are almost 720 exactly the same; for example with H = 1, the spherical harmonic coefficients of Q -(2-H) are within 3% for all n≥1 and the ratio converges rapidly to 1 for large n. The conclusion is that filtering the anomaly by and the multiplying the mean by the above factor is a practical method of fractionally integrating a function on the sphere.