Milestone 5 - Planetary Heat Transfer

  1. Transport in the Atmosphere and Oceans
  2. Modeling Turbulent and Polewards Heat Transport with Diffusion
  3. Diffusion Operator in Spherical Coordinates
  4. Choice of the Diffusion Coefficient in the 2D EBM

Transport in the Atmosphere and Oceans

As a next improvement of our EBM, we want to connect the temperatures of the each individual grid cell to form a globally coupled system that models a smooth surface temperature field of the whole Earth. If we consider a fully coupled system where all grid cells are connected (and not independently solved as tested in milestone 4), heat will get transported from warm regions to colder regions, i.e., heat will be transported polewards. The distribution of the solar forcing term clearly shows that there is more net solar radiation in tropical latitudes than in the polar regions. Thus, some heat transfer polewards is expected. On Earth, the atmosphere and the oceans include mechanisms and processes that cause such a poleward heat transport. For the interested reader, we refer for instance to the following introductory videos that discuss some of the important processes.

  • Global atmospheric circulation on Earth

  • Effect of Coreolis force in large-scale vortices



An important property of moving fluid flows is turbulent motion. Most flows in nature (and engineering) are turbulent, which means multiscale in space and time. From very small eddies in the range of milimeters up to large vortical structures of size hundreds or thousands of meters. Turbulent flows are highly sensitive to small scale disturbances (this is known as the butterfly effect), almost chaotic and thus hard to predict. One of the very important mechanisms of turbulent flows is increased mixing efficiency of cold and warm fluids, i.e., an enhanced heat transfer and heat distribution. This formally requires the numerical resolution (very small grid cells) of all small and large eddies. Small grid cells of course make the simulations very computationally intense and unfortunately for realistic flow scenarios too expensive, even on today's largest super computers. Hence, the full resolution of turbulent fluids is an unsolved problem in fluid mechanics. The current strategy in computational fluid dynamics (and also atmospheric science) to deal with this problem is to introduce additional parametrizations (also often coined turbulence models) that account for the effects of small eddies that cannot be resolved by the grid, so-called subgrid effects. Hence, all modern ESM include a full atmosphere simulation model, which itself includes parametrization for subgrid scale turbulence.

Bad news: Unfortunately, a full atmospheric model is out of the scope of the course.
Good news: It is possible to model the increased mixing of the turbulent fluid and the general convection transport of heat towards the pole by a simplified (but highly parametrized) process, namely by a diffusion process.

Modeling Turbulent and Polewards Heat Transport with Diffusion

In his pioneering work, Budyko

Budyko, M. I. (1969). The effect of solar radiation variations on the climate of the Earth. tellus, 21(5), 611-619.

considered a simple 1D (latitude dependent) EBM that included a simple term that models heat transfer accross latitude zones. From satellite data one can reconstruct that the main heat transfer depends on the latitude and is polewards. The next figure, for instance, shows the heat transfer in the northern hemisphere. It is also worth noting that the heat transfer in the atmosphere is different compared to the ocean.

  • Heat transfer in the northern hemisphere. Generated with data from "Maslin (2013). 'Climate: A Very Short Introduction'"

The idea of Budyko was to include an algebraic term that is proportional to

∼D [T(θ)−Tavg], \sim D\,[T(\theta) - T_{avg}],

where T(θ)T(\theta) is the temperature at latitude θ\theta, TavgT_{avg} is the mean global surface temperature, and DD is a parameter that can be adjusted. The next figure shows a sketch of the model with energy fluxes across the latitude zones (rings), where the coefficient DD is denoted with CC in this figure.

Sellers,

Sellers, W. D. (1969). A global climatic model based on the energy balance of the earth-atmosphere system. Journal of Applied Meteorology and Climatology, 8(3), 392-400.

proposed a model for the heat transport that is based on a diffusion operator. Diffusion or heat conduction is a model where the heat flux is proportional to the gradient of the temperature

∼D ∇⃗T, \sim D\,\vec{\nabla} T,

and the diffusion operator is of the form

∼∇⃗⋅(D ∇⃗T), \sim \vec{\nabla}\cdot (D\,\vec{\nabla} T),

where DD is now the so-called diffusion coefficient. The diffusion coefficient needs to include the effect of oceanic heat transport, turbulent mixing, the sensible heat transport and the latent heat transport. Sensible heat transfer occurs when objects with different temperatures are in direct contact (e.g. surface of the land with the atmosphere) and exchange heat. Latent heat cannot be measured directly, but is an important form of heat transfer that is relevant when a material changes its state (from liquid to gas or from ice to liquid).

Remark: It is clear that all these complex processes cannot be phenomenologically modeled by a single parameter DD. Coming back to our discussion on the different modeling strategies, we already had models based on first principles. These models are (supposed to be) universally valid and based on fundamental laws of nature with no tunable parameters. A next step was to approximate processes that can be measured with simplified models. These empirical models are tuned with parameters to approximate the behavior of the processes observed in measurements. The empircal models are still grounded in reality as they are calibrated with measurements, but include some for of parameter tuning. As an example, we had the radiation model of Budyko in milestone 2.

Due to the complexity of the heat transport and the many mechanisms and processes involved, we do not consider first principle modeling, nor empirical models for the heat transfer in our EBM. Instead, we resort now to another strategy in modeling. We add a simple replacement model, which does not directly mimic the processes, but mimics the effect of the processes on the solution of the EBM. To be precise, we aim to model the effect of the poleward heat transfer on the resulting EBM temperature. This simple replacement model has again parameters, which cannot be tuned to measurements and observations of processes. Instead, the idea is to tune the parameters of the model such, that the resulting solution of the model (the surface temperature) behaves as realistic as possible.

It is important to understand that this type of modeling is driven by data of the desired solution (temperature distribution), i.e., we fit the model to existing solutions (for instance from observations) as good as possible. So we do not tune the model of the process, but tune the outcome of the whole model. It is thus clear that such type of modeling is highly heuristic and should be used with extreme caution, as its predictive power is hard to gauge, i.e., expert knowledge and validation is necessary to give such an approach credibility.

Diffusion Operator in Spherical Coordinates

As we solve our EBM on the surface of Earth, i.e., in spherical coordinates, we need to define the gradient operators and the divergence operators in our natural coordinate system.

The gradient and divergence operators defined in Cartesian coordiates read as,

∇⃗T:=∂T∂xx^+∂T∂yy^+∂T∂zz^,∇⃗⋅F⃗:=∂Fx∂x+∂Fy∂y+∂Fz∂z,\begin{aligned} \vec{\nabla} T &:= \frac{\partial { T}}{\partial { x} } \hat{x} + \frac{\partial { T}}{\partial { y} } \hat{y} + \frac{\partial { T}}{\partial { z} } \hat{z}, \\ \vec{\nabla} \cdot \vec{F} &:= \frac{\partial { F_x}}{\partial { x} } + \frac{\partial { F_y}}{\partial { y} } + \frac{\partial { F_z}}{\partial { z} }, \end{aligned}

where TT is a scalar field, x^\hat{x}, y^\hat{y}, and z^\hat{z} are the unit vectors in the xx, yy, and zz directions, respectively, and F⃗=Fxx^+Fyy^+Fzz^=(Fx,Fy,Fz)\vec{F}=F_x \hat{x} + F_y \hat{y} + F_z \hat{z} = (F_x, F_y, F_z) is a vector field.

We can transform the gradient and divergence operators to the spherical coordinate system by computing the relation between the unit vectors and partial derivatives of both coordinate systems.

For instance, for the unit vectors we have

x^=sin⁡θ~cos⁡φr^+cos⁡θ~cos⁡φθ~^−sin⁡φφ^y^=sin⁡θ~sin⁡φr^+cos⁡θ~sin⁡φθ~^+cos⁡φφ^z^=cos⁡θ~r^−sin⁡θ~θ~^,\begin{aligned} \hat{x} &= \sin \tilde \theta \cos \varphi \hat{r} + \cos \tilde \theta \cos \varphi \hat{\tilde \theta} - \sin \varphi \hat{\varphi} \\ \hat{y} &= \sin \tilde \theta \sin \varphi \hat{r} + \cos \tilde \theta \sin \varphi \hat{\tilde \theta} + \cos \varphi \hat{\varphi} \\ \hat{z} &= \cos \tilde \theta \hat{r}- \sin \tilde \theta \hat{\tilde \theta}, \end{aligned}

and for the partial derivatives we have

∂∂x=sin⁡θ~cos⁡φ∂∂r+cos⁡θ~cos⁡φr∂∂θ~−sin⁡φrsin⁡θ~∂∂φ∂∂y=sin⁡θ~sin⁡φ∂∂r+cos⁡θ~sin⁡φr∂∂θ~+cos⁡φrsin⁡θ~∂∂φ∂∂z=cos⁡θ~∂∂r−sin⁡θ~r∂∂θ~.\begin{aligned} \frac{\partial { }}{\partial { x} } &= \sin \tilde \theta \cos \varphi \frac{\partial { }}{\partial { r} } + \frac{\cos \tilde \theta \cos \varphi}{r} \frac{\partial { }}{\partial { \tilde \theta} } - \frac{\sin \varphi}{r \sin \tilde \theta} \frac{\partial { }}{\partial { \varphi} } \\ \frac{\partial { }}{\partial { y} } &= \sin \tilde \theta \sin \varphi \frac{\partial { }}{\partial { r} } + \frac{\cos \tilde \theta \sin \varphi}{r} \frac{\partial { }}{\partial { \tilde \theta} } + \frac{\cos \varphi}{r \sin \tilde \theta} \frac{\partial { }}{\partial { \varphi} } \\ \frac{\partial { }}{\partial { z} } &= \cos \tilde \theta \frac{\partial { }}{\partial { r} } - \frac{\sin \tilde \theta}{r} \frac{\partial { }}{\partial { \tilde \theta} }. \end{aligned}

Inserting (5) and (6) into (4), we obtain the gradient and divergence operators in spherical coordinates,

∇⃗T=∂T∂rr^+1rsin⁡(θ~)∂T∂φφ^+1r∂T∂θ~θ~^∇⃗⋅F⃗=1r2∂∂r(Frr2)+1rsin⁡(θ~)∂Fφ∂φ+1rsin⁡(θ~)∂∂θ~(Fθ~sin⁡(θ~))\begin{aligned} \vec{\nabla} T &= \frac{\partial T}{\partial r} \hat{r} + \frac{1}{r\sin(\tilde \theta)}\frac{\partial T}{\partial \varphi} \hat{\varphi} + \frac{1}{r}\frac{\partial T}{\partial \tilde \theta} \hat{\tilde \theta} \\ \vec{\nabla} \cdot \vec{F} &= \frac{1}{r^{2}}\frac{\partial { }}{\partial { r} }({F}_r r^{2}) + \frac{1}{r\sin(\tilde \theta)} \frac{\partial {F}_{\varphi}}{\partial \varphi} + \frac{1}{r\sin(\tilde \theta)}\frac{\partial}{\partial \tilde \theta} ({F}_{\tilde \theta}\sin(\tilde \theta)) \end{aligned}

The gradient operator is also commonly written in the vector form (r,φ,θ~)(r,\varphi,\tilde \theta) as

∇⃗T=(∂T∂r,1rsin⁡(θ~)∂T∂φ,1r∂T∂θ~),\begin{aligned} \vec{\nabla} T = \left(\frac{\partial T}{\partial r}, \frac{1}{r\sin(\tilde \theta)}\frac{\partial T}{\partial \varphi}, \frac{1}{r}\frac{\partial T}{\partial \tilde \theta}\right), \end{aligned}

We obtain the diffusion operator in spherical coordinates by combining the expressions for the gradient and divergence operators to obtain

∇⃗⋅(D∇⃗T)=1r2∂∂r(Dr2∂T∂r)+csc⁡2(θ~)r2∂∂φ(D∂T∂φ)+csc⁡(θ~)r2∂∂θ~(Dsin⁡(θ~)∂T∂θ~)\begin{aligned} \vec{\nabla} \cdot (D\vec{\nabla} T) = \frac{1}{r^2} \frac{\partial { }}{\partial { r} } \left( D r^2 \frac{\partial { T}}{\partial { r} } \right) + \frac{\csc^{2}(\tilde \theta)}{r^2} \frac{\partial}{\partial \varphi}\Bigl(D\frac{\partial T}{\partial \varphi}\Bigr) + \frac{\csc (\tilde \theta)}{r^2} \frac{\partial}{\partial \tilde \theta}\left(D \sin (\tilde \theta)\frac{\partial T}{\partial \tilde \theta}\right) \end{aligned}

Note that these derivations are in full three-dimensional space. However, in our EBM, we consider only the surface of the Earth, i.e. a 2D2D approximation of the surface temperature with the choice of the radius coordinate r=REr = R_E. Hence, heat transfer in vertical direction (along the radial coordinate rr) is neglected and only the variations in latitude and longitude directions are considered. Furthermore, we absorb the scaling with the Earth radius RER_E into the diffusion coefficient and simplify the expression to obtain

∇⃗⋅(D∇⃗T)=csc⁡2(θ~)∂∂φ(D~∂T∂φ)⏟Term 1+∂∂θ~(D~∂T∂θ~)⏟Term 2+cot⁡(θ~)D~∂T∂θ~.⏟Term 3\begin{aligned} \vec{\nabla} \cdot (D\vec{\nabla} T) = \underbrace{ \csc^{2}(\tilde \theta) \frac{\partial}{\partial \varphi}\Bigl(\widetilde D\frac{\partial T}{\partial \varphi}\Bigr) }_{\text{Term}~1} + \underbrace{ \frac{\partial}{\partial \tilde \theta}\Bigl(\widetilde D\frac{\partial T}{\partial \tilde \theta}\Bigr) }_{\text{Term}~2} + \underbrace{ \cot(\tilde \theta)\widetilde D\frac{\partial T}{\partial \tilde \theta}. }_{\text{Term}~3} \end{aligned}
Remark: The diffusion coefficient that we use in our 2D2D EBM model is scaled with the Earth radius: D~:=D/RE2\widetilde D := D / R_E^2.

Choice of the Diffusion Coefficient in the 2D EBM

Our final form of the 2D2D EBM with heat diffusion reads as

C(x)∂T∂t+A(CO2)+BT−∇⃗⋅(D∇⃗T)=Ssol(x,t). C(x) \frac{\partial { T}}{\partial { t} } + A(CO_2) + B T - \vec{\nabla} \cdot (D\vec{\nabla} T) = S_{sol}(x,t).

respectively, in spherical coordinates

C(θ~,φ)∂T∂t+A(CO2)+BT−[csc⁡2(θ~)∂∂φ(D~∂T∂φ)+∂∂θ~(D~∂T∂θ~)+cot⁡(θ~)D~∂T∂θ~]=Ssol(θ~,φ,t).\begin{aligned} C(\tilde \theta,\varphi) \frac{\partial { T}}{\partial { t} } &+ A(CO_2) + B T \\ &- \left[ \csc^{2}(\tilde \theta) \frac{\partial}{\partial \varphi}\Bigl(\widetilde D\frac{\partial T}{\partial \varphi}\Bigr) + \frac{\partial}{\partial \tilde \theta}\Bigl(\widetilde D\frac{\partial T}{\partial \tilde \theta}\Bigr) + \cot(\tilde \theta)\widetilde D\frac{\partial T}{\partial \tilde \theta}\right] = S_{sol}(\tilde \theta,\varphi,t). \end{aligned}
Remark: An analysis of the physical dimension of the diffusion coefficient shows that DD has the SI units [WK]\left[\frac{W}{K}\right] (Watts per Kelvin), and D~\widetilde D has the SI units [Wm2K]\left[\frac{W}{m^2 K}\right] (Watts per Kelvin per square meter).

As discussed above, the heat diffusion term needs to model the complex mechanism of poleward transport of heat. There are many different choices of DD available in literature, from simple global constant values to more complex approximations. Here, again, we follow the paper by Zhuang et al.

Zhuang, K., North, G. R., & Stevens, M. J. (2017). A NetCDF version of the two-dimensional energy balance model based on the full multigrid algorithm. SoftwareX, 6, 198-202.

We distinguish between oceanic heat transport and heat transport over land and snow/ice covered areas. We also account for a difference of the heat transfer in the northern and southern hemisphere due to the asymmetric distribution and sizes of the land masses. The following figure shows a sketch of the energy fluxes due to heat transfer accounted in our model (here, the diffusion coefficients are denoted with KK instead of D~\widetilde D)

Thus, for the actual values of the diffusion coefficients we differentiate between grid cells that represent ocean, and all other types. Furthermore, for non oceanic grid cells, we differentiate if we are in the northern or southern hemisphere.

The values for oceanic grid cells in physical unit [W/m2/K][W/m^2/K] are

D~=D~ocean,poles+(D~ocean,equ−D~ocean,poles) sin5(θ~), \widetilde D = \widetilde D_{\text{ocean,poles}} + (\widetilde D_{\text{ocean,equ}} - \widetilde D_{\text{ocean,poles}})\,sin^5(\tilde \theta),

with D~ocean,poles=0.4\widetilde D_{\text{ocean,poles}} = 0.4 and D~ocean,equ=0.65\widetilde D_{\text{ocean,equ}} = 0.65.

The values for non oceanic grid cells in the northern hemisphere in physical units [W/m2/K][W/m^2/K] are

D~=D~NP+(D~equ−D~NP) sin5(θ~), \widetilde D = \widetilde D_{\text{NP}} + (\widetilde D_{\text{equ}} - \widetilde D_{\text{NP}})\,sin^5(\tilde \theta),

with D~NP=0.28\widetilde D_{\text{NP}} = 0.28 and D~equ=0.65\widetilde D_{\text{equ}} = 0.65.

The values for non oceanic grid cells in the southern hemisphere in physical units [W/m2/K][W/m^2/K] are

D~=D~SP+(D~equ−D~SP) sin5(θ~). \widetilde D = \widetilde D_{\text{SP}} + (\widetilde D_{\text{equ}} - \widetilde D_{\text{SP}})\,sin^5(\tilde \theta).

with D~SP=0.20\widetilde D_{\text{SP}} = 0.20 and D~equ=0.65\widetilde D_{\text{equ}} = 0.65.

Remark: We stress again, that the choice of the diffusion coefficients is kind of arbitrary and highly tuned to fit the resulting temperature fields as good as possible. Thus, the choice of a ramp up function between the values at the equator and the values at the poles proportional to sin5(θ~)sin^5(\tilde \theta) has no deeper meaning but is a modeling/tuning choice. The resulting distribution of the diffusion coefficients over ocean and non oceanic grid cells is shown in the following figure

Note that this plot is showing the latitude rather than the colatitude.

Created by Gregor Gassner and Andrés Rueda-Ramírez with contributions by Simone Chiocchetti, Daniel Bach, Sophia Horak, Philipp Baasch, Benjamin Bolm, Erik Faulhaber, and Luca Sommer. Last modified: April 02, 2026. Website built with Franklin.jl and the Julia programming language.