A simple thermal model of the planet’s surface

1 Introduction

The Earth absorbs energy from the Sun. Energy accumulates beneath its surface. It releases energy into the cosmos through various heat transfer mechanisms that operate in parallel: conduction, evaporation, convection, and radiation. The contribution of each of these mechanisms to the total energy released is poorly understood. All that is known is that the hotter the Earth’s surface, the more energy is released.

It is worth remembering that the Sun is by far the Earth’s main source of energy. Solar radiation is considered to represent 99.97% of the Earth’s energy budget. See here and here.

2 Differential equation of the model

To calculate the energy balance of this diagram, simply note that the change in accumulated energy is equal to the energy absorbed minus the energy dissipated.

Let us consider any point on the Earth’s surface. Assume that at time t, the surface temperature is equal to T, and that it has the value T+dT at time t+dt.

Ignoring energy exchanges with neighboring regions, for a small surface area around this point, we can consider the change in accumulated energy to be proportional to the change in temperature dT. Let E_{abs} and E_{ev} be the energies absorbed and dissipated during the time interval dt, and we have:

$$\begin{array}{}\text{(1)}& CdT=E_{abs}-E_{ev}\end{array}$$

If $P_{abs}$ is the absorbed power and $P_{ev}$ the evacuated power,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{E}_{abs}& =P_{abs}dt\\ E_{ev}& =P_{ev}dt\end{array}\end{array}$$

Eq. 1 becomes

$$\begin{array}{}\text{(3)}& \frac{dT}{dt}=\frac{1}{C}(P_{abs}-P_{ev})\end{array}$$

Using Stockwell notation

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{P}_{abs}& =S\\ P_{ev}& =F\end{array}\end{array}$$

We get

$$\begin{array}{}\text{(5)}& \frac{dT(t)}{dt}=\frac{1}{C}(S(t)-F(T))\end{array}$$

If the Earth is subjected to constant radiation, a thermal equilibrium is eventually established. At this moment, the power evacuated is equal to the power absorbed, and the temperature remains constant:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\frac{dT(t)}{dt}& =0\\ S& =F\end{array}\end{array}$$

Not much more can be learned from this equation. $S(t)$ depends directly on the solar radiation at the top of the atmosphere, to which a poorly known coefficient must be applied to account for the passage through the atmosphere and reflection at the surface. The $F(T)$ law is very poorly understood. It contains a $T^4$ term corresponding to the power removed by thermal radiation, and one or more lower power terms for other modes of heat removal.

On the other hand, if we focus on perturbations relative to an equilibrium situation, we can learn a lot from them.

Let us therefore assume that at time $t=0$, $T=T_0$, $S=S_0$ and $F=F_0=S_0$. Let us linearize $F(T)$ by a first-order Taylor expansion. If we set $\frac{1}{\lambda }=F^{\prime }(T_0)$ $$\begin{array}{}\text{(7)}& F(T)=F_0+\frac{T-T_0}{\lambda }\end{array}$$

Let us again set $T_a=T-T_0$ (the temperature anomaly with respect to $T_0$), $S_a=S-S_0$ (the radiation anomaly with respect to $S_0$) and noting that $d{T}_a=dT$ and that $F_0=S_0$ , Eq. 7 is written

$$\begin{array}{}\text{(8)}& \frac{d{T}_a(t)}{dt}=\frac{1}{C}(S_a(t)-\frac{T_a(t)}{\lambda })\end{array}$$

Setting $k=\lambda C$, the solution to this equation is given by
$$\begin{array}{}\text{(9)}& T_a(t)=\frac{\lambda }{k}{e}^{-t/k}{\int }_0^t{e}^{u/k}{S}_a(u)du\end{array}$$

It is mathematically proven that this equation is convergent to an asymptotic solution when the coefficient k is positive.

3 Analytic Solutions

3.1 Step Function

If the system is at equilibrium at time $t=0$ at temperature $T_0$ and $S_a(t)=\mathrm{\Delta }{S}_0$ for all $t>0$ , then

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathrm{\Delta }T(t)& =T(t)-T_0\\ & =\frac{\lambda }{k}{e}^{-t/k}k\mathrm{\Delta }{S}_0(e^{t/k}-1)\\ & =\lambda \mathrm{\Delta }{S}_0(1-e^{-t/k})\end{array}\end{array}$$

The temperature variation evolves asymptotically towards the value $\lambda \mathrm{\Delta }{S}_0$. See Fig. 1.

Warning: Detected the use of `TeX()`, but mathjax has not been specified. Try
running `config(.Last.value, mathjax = 'cdn')`

Fig. 1: Eq. 10 with $\lambda =\mathrm{\Delta }{S}_0=1$, for different values of k

3.2 General function

A general function $S_a(t)$ can be discretized into several staircases. If we assume it to be constant in the intervals $i,i+1$, using Eq. 10, we obtain

$$\begin{gathered} \begin{array}{}\text{(11)}& T(t)-T_0=\lambda \sum _i\mathrm{\Delta }{S}_i(1-e^{-(t-i)/k}) \\ \mathit{with} \\ \mathrm{\Delta }{S}_i=S_{i+1}-S_i=S_{a_{i+1}}-S_{a_i}\end{array} \end{gathered}$$

If we approximate $(1-e^{-t/k})$ by a linear-step function equal to $t/\tau$ for $t$ between 0 and $\tau$, and equal to 1 for t greater than $\tau$ (see Fig. 2), Eq. 11 can be greatly simplified. Indeed, the formula can then be written in finite difference form:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}T(t)-T_0& =\lambda (T_1+T_2)\\ T_1& =\sum _{i=0,1,...,t-\tau }\mathrm{\Delta }{S}_i\\ T_2& =\sum _{i=0,1,...,\tau }\mathrm{\Delta }{S}_{t-\tau +i}\frac{\tau -i}{\tau }\end{array}\end{array}$$

Noting

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\overline{S}}_{t-\tau }^t& =theaverageofSbetweentimestandt-\tau \\ {\overline{S}}_{a_{t-\tau }}^t& =theaverageof{S}_{a}betweentimestandt-\tau \end{array}\end{array}$$

By expanding the sums and grouping the terms, we easily obtain

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{T}_1& =S_{t-\tau }-S_0\\ T_2& ={\overline{S}}_{t-\tau }^t-S_{t-\tau }\end{array}\end{array}$$

It comes finally

$$\begin{array}{}\text{(15)}& T(t)-T_0=\lambda ({\overline{S}}_{t-\tau }^t-S_0)=\lambda {\overline{S}}_{a_{t-\tau }}^t\end{array}$$

Fig. 2: Example of approximation of an exponential function by a linear function at the start and constant thereafter.

3.3 Sinusoidal variation

Let us suppose that $S_a(t)$ is a sinusoid of amplitude $\mathrm{\Delta }{S}_a$ and period $per$:

$$\begin{array}{}\text{(16)}& S_a(t)=\mathrm{\Delta }{S}_{a}sin(2\pi \frac{t}{per})\end{array}$$

In this case the solution is given by

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }T(t)=T(t)-T_0=e^{-t/k}Asin(\varphi )+Asin(\omega t-\varphi )\end{array}$$

In that last formula,

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\omega & =\frac{2\pi }{per}\\ \varphi & =atan(k\omega )\\ A& =\frac{1}{\sqrt{1+(k\omega {)}^2}}\lambda \mathrm{\Delta }{S}_a\end{array}\end{array}$$

The first term in Eq. 17 is a transient term. It tends to 0 as $t$ tends to infinity, and $\mathrm{\Delta }T(t)$ tends to a sinusoid out of phase by an angle $\varphi$ with respect to the solar radiation.

For slow oscillations, $k\omega$ tends to zero, the amplitude of the temperature variation tends to $\lambda \mathrm{\Delta }{S}_a$ and $\varphi$ tends to zero.

For fast oscillations, $k\omega$ tends to infinity, the amplitude of the temperature variation tends to zero and $\varphi$ tends to $\pi /2$.

See Fig. 3.

Warning: Detected the use of `TeX()`, but mathjax has not been specified. Try
running `config(.Last.value, mathjax = 'cdn')`

Fig. 3: Temperature anomaly for a sinusoidal solar anomaly (Eq. 17, without the transient term), with $\lambda =\mathrm{\Delta }{S}_a=1,andk=5$, for different values of the period $per$. On the abscissa, the time and on the ordinate the temperature anomaly, without particular units.
It should be noted that the amplitude of the temperature anomaly is very small when the period of the input signal is short. It is thermal inertia that is at work.

3.4 Out-of-phase sinusoidal variation

If $\mathrm{\Delta }S(t)$ is out of phase by an angle $\alpha$:

$$\begin{array}{}\text{(19)}& \mathrm{\Delta }S(t)=\mathrm{\Delta }{S}_{a}sin(2\pi \frac{t}{per}-\alpha )\end{array}$$

In this case the solution of Eq. 9 becomes

$$\begin{array}{}\text{(20)}& \mathrm{\Delta }T(t)=T(t)-T0=e^{-t/k}Asin(\varphi +\alpha )+Asin(\omega t-(\varphi +\alpha ))\end{array}$$

Les paramètres de cette formule sont également calculés selon l’Eq. 18.

3.5 Any periodic function

If $S_a(t)$ is a periodic function, it can be decomposed into a Fourier series:

$$\begin{array}{}\text{(21)}& S_a(t)=\mathrm{\Delta }{S}_0+\sum _{i=1..n}\mathrm{\Delta }{S}_{i}sin({\omega }_{i}t+{\beta }_i)\end{array}$$

The term $\mathrm{\Delta }{S}_0$ corresponds to the 0th harmonic of the decomposition and is equal to the average of $S_a(t)$ over its period. The solution will be obtained by applying Eq. 10 to it and Eq. 19 to the other terms.

3.5.1 Example: incident solar radiation at the spring or autumn equinox

At temperate latitudes, within a scale factor, the radiation is very well modeled by a half sinusoid for the 12 hours following sunrise, and a zero value for the following 12 hours. The mean of this function is equal to $1/\pi$. The anomaly of this radiation can be expressed as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rlr}{S}_a(t)& =\mathrm{\Delta }{S}_0(sin(\omega t)-1/\pi )& 0<=t<=\frac{\pi }{\omega }\\ & =-\mathrm{\Delta }{S}_0/\pi & \frac{\pi }{\omega }<=t<=\frac{2\pi }{\omega }\end{array}\end{array}$$

For $\mathrm{\Delta }{S}_0=1$ , the Fourier transform of this function has an analytical solution given by the following formula:

$$\begin{array}{}\text{(23)}& fft(S_a(t))=\frac{1}{2}sin(\omega t)+\frac{2}{\pi }\sum _{n=2,4,6,....}\frac{1}{n^2-1}sin(n\omega t-\frac{\pi }{2})\end{array}$$

It will be enough to apply the previous formulas to each term. See Fig. 4.

Warning: Detected the use of `TeX()`, but mathjax has not been specified. Try
running `config(.Last.value, mathjax = 'cdn')`

Fig. 4: It takes the theoretical solar radiation on an equinox day. This radiation is approximated by a half-sinusoid for 12 hours after sunrise, followed by a zero value during the night. It is represented by the curve Sa.
The temperature anomaly was calculated by applying the formulas of Eq. 18 and Eq. 20 to each term of Eq. 23. The calculations were carried out considering $\lambda =1$, for some values ​​of $k$. On the abscissa, the time expressed in hours since sunrise. On the ordinate, the solar radiation anomaly and the temperature anomalies (without particular units).

4 Analogies to the thermal model

4.1 Draining a container supplied with water

Consider a container of constant cross-section supplied with water through a faucet and equipped with a drain hole at its base. The water balance in the container is easy to establish.

Noting:

$S$	The recipient sectiont
$s$	The drain hole section
$h$	The height of water above the drain hole
$D_{in}$	The feed rate
$D_{out}$	The evacuated flow

On a:
$$\begin{array}{}\text{(24)}& \frac{dh}{dt}=\frac{1}{S}(D_{in}-D_{out})\end{array}$$

This is analogous to Eq. 3. Just consider that a volume of water is energy, a flow rate of water is power, and the height of water is temperature.

Using Torricelli’s formula, we easily find that

$$\begin{array}{}\text{(25)}& D_{out}=s\sqrt{2gh}\end{array}$$

At equilibrium, for a constant feed rate, we obtain

$$\begin{array}{}\text{(26)}& h=\frac{1}{2g}(D_{in}/s{)}^2\end{array}$$

Finally, Eq. 24 is written

$$\begin{array}{}\text{(27)}& \frac{dh}{dt}=\frac{1}{S}(D_{in}-s\sqrt{2gh})\end{array}$$

As before, this equation can be linearized with respect to an equilibrium situation. Suppose that at time $t=0$, $h=h_0$, $D_{in}=D_{i{n}_0}$ and $D_{out}=D_{ou{t}_0}=s\sqrt{2g{h}_0},D_{ou{t}_0}^{\prime }=s\sqrt{g/2h_0}$ .

Noting $h^a=h-h_0,D_{in}^a=D_{in}-D_{i{n}_0}$, it finally comes

$$\begin{array}{}\text{(28)}& \frac{d{h}^a}{dt}=\frac{1}{S}(D_{in}^a-h^{a}s\sqrt{g/2h_0})\end{array}$$

This will allow us to judge to what extent the linearization approximation deviates from the exact solution.

4.1.1 A concrete example for a sinusoidal power supply flow

4.1.1.1 Exact solution

Eq. 24 was numerically integrated for a container with a cross-section $S=100c{m}^2$, equipped with a hole with a cross-section $s=0.5c{m}^2$ fed by a variable flow rate $D_{in}=D_0(1+ksin(2\pi t/per))$ with $D_0=3l/min$. The initial water height $h_0$ was calculated by Eq. 26 with $D_{in}=D_0$.

Fig. 5: Solution of Eq. 24 by numerical integration for a concrete case, for $k=0.5$ and some values ​​of $per$. We can observe the convergence towards a sinusoid which oscillates around the average value of the input flow. The amplitude and the phase shift of the sinusoid are analogous to what we observe in Fig. 3

4.1.1.2 Approximate solution by linearization

Fig. 6: Solution of Eq. 28 by numerical integration for a concrete case, for $k=0.5$ and some values ​​of $per$. The differences with the Fig. 5 are almost imperceptible.

Fig. 7: Approximation errors by linearization are negligible.

4.1.2 Response to an impulse

Suppose the same container as before has reached equilibrium at a feed rate $D_0=3l/min$ , and that for a brief period this flow rate is doubled.

Fig. 8: Exact solution of Eq. 28 by numerical integration, for a doubling of the flow rate over 5 seconds, starting from an equilibrium situation. The level increases rapidly during the doubling before decreasing much more slowly towards the initial equilibrium.

A similar phenomenon is probably at work in the Dansgaard-Oeschger events during which we observe an increase in temperature of 5 to 8°C over a few decades followed by a cooling period extending over several hundred years.

4.2 Equation of motion of a flywheel

Indoor bikes are often equipped with a flywheel that dissipates the power produced by the cyclist through eddy currents. The differential equation for such a system subject to a motor torque and a resistive torque is easy to formulate based on the principle of conservation of energy.

Noting:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\theta & =Angularpositionoftheflywheel\\ \omega =\dot{\theta }& =Angularspeedoftheflywheel\\ \dot{\omega }=\ddot{\theta }& =Angularaccelerationoftheflywheel\\ J& =Rotationalinertiaoftheflywheel\\ E_v=J{\omega }^2/2& =Kineticenergyoftheflywheel\\ C_m& =Drivingtorque\\ C_r& =Resistanttorque\\ P_m=C_m\omega & =Drivingpower\\ P_r=C_r\omega & =Resistantpower\\ E_m=P_{m}dt& =Energyinjectedintotheflywheelduringtheintervaldt\\ E_r=P_{r}dt& =Energylostbytheflywheelduringtheintervaldt\end{array}\end{array}$$

The energy balance of the flywheel is written:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}d{E}_v& =E_m-E_r\\ J\omega d\omega & =P_{m}dt-P_{r}dt\\ J\omega \dot{\omega }& =C_m\omega -C_r\omega \\ \\ \dot{\omega }=\ddot{\theta }& =\frac{1}{J}(C_m-C_r)\end{array}\end{array}$$

This equation is equivalent to Eq. 5. It allows us to better understand the meaning that must be given to thermal inertia.

5 Determining model parameters based on observations

Several temperature series are available. They are global or local. Some series are published in grid form, covering all or part of the Earth’s surface. If the temperatures are absolute, it will be easy to derive an anomaly by difference with a mean.

Absorbed radiation is poorly understood. At some stations, incident radiation has been measured for many years, but not reflected radiation.

In the absence of local measurements of incident radiation, it will be necessary to rely on total solar irradiance (TSI). This has only been measured by satellite for a few decades, but historical reconstructions of the TSI exist that go back much further. When the calculation is performed at a specific point, it will be possible to convert the TSI to the theoretical local incident radiation using astronomical calculations. Specialized libraries are available that greatly facilitate these calculations.

Eq. 8, Eq. 9, and Eq. 15 can be used in all cases, whether solar activity is periodic or not.

Eq. 8 requires numerical integration, which can be quite cumbersome.

Eq. 9 is not recommended because multiplying very small numbers by very large numbers will generally cause related numerical instabilities.

Eq. 15 is fairly straightforward to use, but requires knowledge of the solar activity history over a sufficiently long period before the start of the temperature series.

For periodic solar activity (daily or seasonal), the simplest solution is to decompose the signal into a Fourier series.

In all cases, nonlinear multiple regression will be required to obtain the model parameters.