Warning: Detected the use of `TeX()`, but mathjax has not been specified. Try
running `config(.Last.value, mathjax = 'cdn')`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:
\[ C \ dT = E_{abs} - E_{ev} \tag{1}\]
If \(P_{abs}\) is the absorbed power and \(P_{ev}\) the evacuated power,
\[ \begin{align} E_{abs} & = P_{abs}\ dt \\ E_{ev} & = P_{ev}\ dt \end{align} \tag{2}\]
Eq. 1 becomes
\[ \frac{dT}{dt} = \frac{1}{C}(P_{abs} - P_{ev}) \tag{3}\]
Using Stockwell notation
\[ \begin{align} P_{abs} & = S \\ P_{ev} & = F \end{align} \tag{4}\]
We get
\[ \frac{dT(t)}{dt} = \frac{1}{C}(S(t) - F(T)) \tag{5}\]
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{align} \frac{dT(t)}{dt} & = 0 \\ S & = F \end{align} \tag{6}\]
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'(T_0)\) \[ F(T)=F_0 + \frac{T-T_0}{\lambda} \tag{7}\]
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 \(dT_a=dT\) and that \(F_0=S_0\) , Eq. 7 is written
\[ \frac{dT_a(t)}{dt} = \frac{1}{C}(S_a(t) - \frac{T_a(t)}{\lambda}) \tag{8}\]
Setting \(k=\lambda C\), the solution to this equation is given by
\[
T_a(t) = \frac{\lambda}{k}\ e^{-t/k}\int_{0}^{t} e^{u/k}\ {S_a(u)} \ du
\tag{9}\]
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) = \Delta S_0\) for all \(t>0\) , then
\[ \begin{align} \Delta T(t) & = T(t) - T_0 \\ & = \frac{\lambda}{k} \ e^{-t/k} \ k \ \Delta S_0 \ (e^{t/k} - 1) \\ & = \lambda \ \Delta S_0 \ (1 - e^{-t/k}) \end{align} \tag{10}\]
The temperature variation evolves asymptotically towards the value \(\lambda \Delta S_0\). See Fig. 1.
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{equation} T(t)\ -\ T_0\ =\ \lambda \ \sum _i\ \Delta S_i\ (1\ -\ e^{-(t-i)/k}) \\ \mathit{\ with\ } \\ \Delta S_i=S_{i+1}-S_i = S_{a_{i+1}} - S_{a_i} \end{equation} \tag{11}\]
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{align} T(t) - T_0 &= \lambda(T_1+T_2) \\ T_1\ &=\ \sum _{i=0,1,...,t-\tau}\Delta S_i \\ T_2\ &=\ \sum _{i=0,1,...,\tau}\Delta S_{t-\tau +i}\ \frac{\tau \ -\ i}{\tau } \end{align} \tag{12}\]
Noting
\[ \begin{align} \bar S^{\ t}_{t-\tau } &= the\ average\ of\ S\ between\ times\ t\ and\ t-\tau \\ \bar S^{\ t}_{a_{t-\tau }} &= the\ average\ of\ S_a\ between\ times\ t\ and\ t-\tau \\ \end{align} \tag{13}\]
By expanding the sums and grouping the terms, we easily obtain
\[ \begin{align} T_1\ &=\ S_{t-\tau }\ -\ S_0 \\ T_2\ &=\ \bar S^{\ t}_{t-\tau }\ -\ S_{t-\tau } \end{align} \tag{14}\]
It comes finally
\[ T(t)\ -\ T_0\ =\ \lambda \ (\bar S^{\ t}_{t-\tau }\ -\ S_0) = \lambda\ \bar S^{\ t}_{a_{t-\tau }} \tag{15}\]
3.3 Sinusoidal variation
Let us suppose that \(S_a(t)\) is a sinusoid of amplitude \(\Delta S_a\) and period \(per\):
\[ S_a(t) = \Delta S_a \ sin(2\pi \frac{t}{per}) \tag{16}\]
In this case the solution is given by
\[ \Delta T(t) = T(t) - T_0 = e^{-t/k} \ A \ sin(\phi) + A \ sin(\omega t - \phi) \tag{17}\]
In that last formula,
\[ \begin{align} \omega & = \frac{2\pi}{per} \\ \phi & = atan(k\omega) \\ A & = \frac{1}{\sqrt{1+(k\omega)^2}} \ \lambda \ \Delta S_a \end{align} \tag{18}\]
The first term in Eq. 17 is a transient term. It tends to 0 as \(t\) tends to infinity, and \(\Delta T(t)\) tends to a sinusoid out of phase by an angle \(\phi\) with respect to the solar radiation.
For slow oscillations, \(k\omega\) tends to zero, the amplitude of the temperature variation tends to \(\lambda \Delta S_a\) and \(\phi\) tends to zero.
For fast oscillations, \(k\omega\) tends to infinity, the amplitude of the temperature variation tends to zero and \(\phi\) 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')`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 \(\Delta S(t)\) is out of phase by an angle \(\alpha\):
\[ \Delta S(t) = \Delta S_a \ sin(2\pi \frac{t}{per} - \alpha) \tag{19}\]
In this case the solution of Eq. 9 becomes
\[ \Delta T(t)=T(t) - T0 = e^{-t/k} \ A \ sin(\phi+\alpha) + A \ sin(\omega t - ( \phi+\alpha)) \tag{20}\]
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:
\[ S_a(t) = \Delta S_0 + \sum_{i=1..n} \Delta S_i \ sin(\omega_i t + \beta_i) \tag{21}\]
The term \(\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{align} S_a(t) & = \Delta S_0 \ (sin(\omega t) - 1/\pi) & 0 <= t <= \frac{\pi}{\omega} \\ & = -\Delta S_0/\pi & \frac{\pi}{\omega} <= t <= \frac{2\pi}{\omega} \end{align} \tag{22}\]
For \(\Delta S_0 = 1\) , the Fourier transform of this function has an analytical solution given by the following formula:
\[ 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}) \tag{23}\]
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')`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:
\[
\frac{dh}{dt} = \frac{1}{S}(D_{in} - D_{out})
\tag{24}\]
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
\[ D_{out}=s\sqrt{2gh} \tag{25}\]
At equilibrium, for a constant feed rate, we obtain
\[ h=\frac{1}{2g}(D_{in}/s)^2 \tag{26}\]
Finally, Eq. 24 is written
\[ \frac{dh}{dt} = \frac{1}{S}(D_{in} - s\sqrt{2gh}) \tag{27}\]
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_{in_0}\) and \(D_{out}=D_{out_0}= s\sqrt{2gh_0}, D'_{out_0}=s\sqrt{g/2h_0}\) .
Noting \(h^a=h-h_0, D^a_{in}=D_{in}- D_{in_0}\), it finally comes
\[ \frac{dh^a}{dt} = \frac{1}{S}(D^a_{in} - h^as\sqrt{g/2h_0}) \tag{28}\]
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=100\ cm^2\), equipped with a hole with a cross-section \(s=0.5\ cm^2\) fed by a variable flow rate \(D_{in}=D_0(1+k\ sin(2\pi t/per))\) with \(D_0=3\ l/min\). The initial water height \(h_0\) was calculated by Eq. 26 with \(D_{in}=D_0\).
4.1.1.2 Approximate solution by linearization
4.1.2 Response to an impulse
Suppose the same container as before has reached equilibrium at a feed rate \(D_0=3\ l/min\) , and that for a brief period this flow rate is doubled.
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{align} \theta &= Angular\ position\ of\ the\ flywheel \\ \omega = \dot \theta & = Angular\ speed\ of\ the\ flywheel \\ \dot\omega=\ddot \theta & = Angular\ acceleration\ of\ the\ flywheel \\ J & = Rotational\ inertia\ of\ the\ flywheel \\ E_v = J \omega^2/2 & = Kinetic\ energy\ of\ the\ flywheel \\ C_m & = Driving\ torque \\ C_r & = Resistant\ torque \\ P_m = C_m \omega & = Driving\ power \\ P_r = C_r \omega & = Resistant\ power \\ E_m = P_m dt &= Energy\ injected\ into\ the\ flywheel\ during\ the\ interval\ dt\\ E_r = P_r dt &= Energy\ lost\ by\ the\ flywheel\ during\ the\ interval\ dt\\ \end{align} \tag{29}\]
The energy balance of the flywheel is written:
\[ \begin{align} dE_v &= E_m - E_r \\ J \omega d\omega &= P_m dt - P_rdt \\ J\omega\dot\omega&=C_m\omega - C_r\omega\\ \\ \dot \omega = \ddot \theta &= \frac{1}{J}(C_m - C_r) \\ \end{align} \tag{30}\]
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.