A theoretical model for realistic local climates

TL;DR

This paper introduces a nonlinear theoretical model that predicts local surface climate conditions, accurately reproducing features like seasonal lag and daily temperature asymmetries, adaptable to different planetary parameters.

Contribution

The paper presents a new nonlinear climate model incorporating solar forcing, land-ocean differences, and humidity, with validation against real Earth climate data.

Findings

Model accurately reproduces Earth's climate features

Parameters aligned with Köppen climate zones improve accuracy

Adaptable to planets with different astronomical parameters

Abstract

We write a nonlinear model that predicts the climate (temperature and humidity) on the surface of a small region on Earth, perform numerical investigations using the model, and compare the results to real climate on a variety of regions on Earth. It the parameters are chosen keeping into consideration the climatic K\"oppen zone to which the region belongs, the numerical model accurately reproduces the real climate. The model takes into account the doubly-periodic forcing of the solar radiation (annual and daily), the laws of irradiance, the fact that the Earth has land and oceans with different thermic inertia, and the humidity of the air due to evaporation. This enables us to reproduce remarkable features of Earth's climate such as lag of seasons, lag of noons, and asymmetric evolution of daily temperatures. The model can easily be adapted to planets with non-terrestrial astronomic parameters. We conclude this article with an investigation of an Earth with eccentricity higher than real.

Figures

Tables

Table 1. Table 1: Fundamental astronomical parameters for planet Earth.

$e$	$0.0167$	Eccentricity of Earth’s orbit
$a$	$1.496\times {10}^{11}m$	Average earth-sun distance
$\gamma$	$23.437°$	Earth’s mean obliquity
$\delta$	$-12.8219°$	Angle between solstices and perihelion/aphelion
$D$	$8.616409\times {10}^{4}s$	Period of rotation of the Earth
$Y$	$3.15569\times {10}^{7}s$	Period of revolution of the Earth

Table 2. Table 2: Fundamental physical parameters for the solar irradiance.

Parameter	Value	Description
$\sigma$	$5.670\times {10}^{-8}{\mathrm{Js}}^{-1}{K}^{-4}{m}^{-2}$	Stefan-Boltzmann constant
$R_S$	$6.955\times {10}^{8}m$	Solar radius
$T_S$	$5778K$	Sun superficial temperature

Table 3. Table 3: Reflectance r 𝑟 r , absorbance α 𝛼 \alpha and transmittance τ 𝜏 \tau of solar radiation for various thermodynamic bodies on Earth.

	Reflectance $r$	Absorbance $\alpha$	Transmittance $\tau$
Atmosphere (solar light)	0.23	$0.25$	0.52
Soil	0.2	0.8	0
Desert	0.4	0.6	0
Ocean	0.15	0.85	0
Ice over land	0.85	0.15	0
Ice over oceans	0.62	0.38	0

Table 4. Table 4: Atmospheric absorbance α 0 T superscript subscript 𝛼 0 𝑇 \alpha_{0}^{T} to thermal radiation and emissivity ε 𝜀 \varepsilon for various thermodynamic bodies on Earth.

Parameter	Value	Description
${\alpha }_0^T$	[0.8, 0.95]	Atmosphere absorbance (infrared light)
${\epsilon }_0^d$	0.8	Atmosphere emissivity downwards
${\epsilon }_0^u$	0.45	Atmosphere emissivity upwards
${\epsilon }_1^{soil}$	0.94	Soil emissivity
${\epsilon }_1^{sand}$	0.75	Sand emissivity
${\epsilon }_1^{ice}$	0.85	Ice over land emissivity
${\epsilon }_2^{water}$	0.96	Water emissivity
${\epsilon }_2^{ice}$	0.85	Ice over oceans emissivity

Table 5. Table 5: Heat transfer coefficients between air and land and air and water.

Parameter	Value	Description
$h_{01}$	$[5,40]{\mathrm{Js}}^{-1}{m}^{-2}$	land-air heat transfer coefficient
$h_{02}$	$[5,40]{\mathrm{Js}}^{-1}{m}^{-2}$	water-air heat transfer coefficient

Table 6. Table 6: Geothermal powers.

Parameter	Value	Description
${\eta }_1$	$0.345{\mathrm{Js}}^{-1}{m}^{-2}$	continental geothermal power
${\eta }_2$	$0.802{\mathrm{Js}}^{-1}{m}^{-2}$	oceanic geothermal power

Table 7. Table 7: Evaporation rate, rainfall rate, and latent heat.

Parameter	Value	Description
${\mu }_2$	$5.7\times {10}^{-5}{s}^{-1}$	Evaporation rate from water
${\mu }_1$	$[{10}^{-6},{10}^{-4}]s^{-1}$	Evaporation rates from land
$\nu$	$[1,5]\times {10}^{-5}{s}^{-1}$	Rainfall rate
$\lambda$	$2.26\times {10}^6{\mathrm{Jkg}}^{-1}$	Latent heat of evaporation and condensation

Table 8. Table 8: Heat capacities.

Parameter	Value	Description
$C_0^d$	$3.5\times {10}^5{\mathrm{JK}}^{-1}{m}^{-2}$	dry air thermal capacity
$C_0^v$	$1\times {10}^6{\mathrm{JK}}^{-1}{m}^{-2}$	water vapor thermal capacity
$C_1^{soil}$	$1.0\times {10}^6{\mathrm{JK}}^{-1}{m}^{-2}$	soil thermal capacity
$C_1^{ice}$	$1.0\times {10}^6{\mathrm{JK}}^{-1}{m}^{-2}$	ice thermal capacity
$C_1^{forest}$	$1.7\times {10}^6{\mathrm{JK}}^{-1}{m}^{-2}$	forest thermal capacity
$C_1^{sand}$	$3.2\times {10}^6{\mathrm{JK}}^{-1}{m}^{-2}$	sand thermal capacity
$C_2$	$[1.7,\mathrm{\hspace{0.17em}2.5}]\times {10}^8{\mathrm{JK}}^{-1}{m}^{-2}$	ocean thermal capacity

Equations

$$\dot{T}(t)=(-\lambda T(t)+\mu )+(a\cos (\omega t)+b)$$

$$\frac{\mu +b}{\lambda }+\frac{a}{\sqrt{{\lambda }^2+{\omega }^2}}\cos (\omega (t-\tau )),$$

$$\rho (t)=\frac{a(1-e^2)}{1+e\cos (\vartheta (t))},\dot{\vartheta }(t)=\frac{2\pi }{Y{\sqrt{1-e^2}}^3}(1+e\cos (\vartheta (t)){)}^2,\vartheta (0)=0$$

$$\{\begin{array}{l}{f}_1(t)=e_1(\cos \gamma \cos \delta \cos (\Omega t)-\sin \delta \sin (\Omega t))+\\ +e_2(\cos \gamma \sin \delta \cos (\Omega t)+\cos \delta \sin (\Omega t))-e_3\sin \gamma \cos (\Omega t)\\ f_2(t)=e_1(\sin \delta (-\cos (\Omega t))-\cos \gamma \cos \delta \sin (\Omega t))+\\ +e_2(\cos \delta \cos (\Omega t)-\cos \gamma \sin \delta \sin (\Omega t))+e_3\sin \gamma \sin (\Omega t)\\ f_3(t)=e_2\sin \gamma \sin \delta +e_1\sin \gamma \cos \delta +e_3\cos \gamma .\end{array}$$

$$d(t)=-\cos (\vartheta (t))e_1-\sin (\vartheta (t))e_2.$$

$$n(t)=\cos (\psi )\cos \phi f_1(t)+\sin (\psi )\cos \phi f_2(t)+\sin \phi f_3(t),$$

$$n(t)\cdot d(t)=\cos\!\varphi\sin(\delta-\vartheta(t))\sin(\Omega t+\psi)+\\ -\cos(\delta-\vartheta(t))\big{(}\sin\!\gamma\sin\!\varphi+\cos\!\gamma\cos\!\varphi\cos(\Omega t+\psi)\big{)}.$$

$$I=\sigma T_s^4\frac{R_s^2}{{\rho }^2}.$$

$$W(t)=\max \left\{\sigma T_s^4\frac{R_s^2}{\rho (t{)}^2}n(t)\cdot d(t),0\right\}.$$

$$\{\begin{array}{l}\frac{d{Q}_0^{SR}}{dt}={\alpha }_0(1+p{\tau }_0{r}_1+q{\tau }_0{r}_2)W(t)\\ \frac{d{Q}_1^{SR}}{dt}=p{\tau }_0{\alpha }_{1}W(t)\\ \frac{d{Q}_2^{SR}}{dt}=q{\tau }_0{\alpha }_{2}W(t).\end{array}$$

$$\{\begin{array}{l}\frac{d{Q}_0^{TR}}{dt}=\sigma (p{\alpha }_0^T{\epsilon }_1{T}_1^4+q{\alpha }_0^T{\epsilon }_2{T}_2^4-({\epsilon }_0^d+{\epsilon }_0^u)T_0^4)\\ \frac{d{Q}_1^{TR}}{dt}=p\sigma ({\epsilon }_0^d{T}_0^4-{\epsilon }_1{T}_1^4)\\ \frac{d{Q}_2^{TR}}{dt}=q\sigma ({\epsilon }_0^d{T}_0^4-{\epsilon }_2{T}_2^4).\end{array}$$

$$\frac{dQ}{dt}=h(T_2-T_1),$$

$$\{\begin{array}{l}\frac{d{Q}_0^C}{dt}=p{h}_{01}(T_1(t)-T_0(t))+q{h}_{02}(T_2(t)-T_0(t))\\ \frac{d{Q}_1^C}{dt}=-p{h}_{01}(T_1(t)-T_0(t))\\ \frac{d{Q}_2^C}{dt}=-q{h}_{02}(T_2(t)-T_0(t)),\end{array}$$

$$\frac{d{Q}_1^M}{dt}=p{\eta }_1,\frac{d{Q}_2^M}{dt}=q{\eta }_2,\frac{d{Q}_0^M}{dt}=0.$$

$${\mu }_2=\frac{25+19{\overline{v}}_i}{3600}.$$

$$\nu \overline{U}Y{\rho }_0{\ell }_0={\rho }_2\Pi .$$

$$\begin{array}{r}\frac{dU}{dt}(t)=p{\mu }_1(U_s(T_0(t))-U(t))+q{\mu }_2(U_s(T_0(t))-U(t))-\nu U=\\ =(p{\mu }_1+q{\mu }_2)(U_s(T_0(t))-U)-\nu U\end{array}$$

$$U_s(T)=e^{0.0666T-23.96}$$

$$\frac{d{Q}_1^{LE}}{dt}=-p\lambda \max \left({\rho }_0{\ell }_0\frac{dU}{dt},0\right),\frac{d{Q}_2^{LE}}{dt}=-q\lambda \max \left({\rho }_0{\ell }_0\frac{dU}{dt},0\right),$$

$$\frac{d{Q}_0^{LC}}{dt}=-\lambda \min ({\rho }_0{\ell }_0\frac{dU}{dt},0).$$

$$\frac{d{Q}_0}{dt}=C_0(U)\frac{d{T}_0}{dt},\frac{d{Q}_1}{dt}=p{C}_1\frac{d{T}_1}{dt},\frac{d{Q}_2}{dt}=q{C}_2\frac{d{T}_2}{dt}.$$

$$C_0(U)=c_0^d{\rho }_0^d{\ell }_0+c_0^v{\rho }_0^d{\ell }_{0}U=C_0^d+C_0^{v}U,$$

$$\{\begin{array}{l}(C_0^d+U{C}_0^v)\frac{d{T}_0}{dt}={\alpha }_0(1+p{\tau }_0{r}_1+q{\tau }_0{r}_2)W+\sigma {\alpha }_0^T(p{\epsilon }_1{T}_1^4+q{\epsilon }_2{T}_2^4)+\\ -\sigma ({\epsilon }_0^d+{\epsilon }_0^u)T_0^4)+p{h}_{01}(T_1-T_0)+q{h}_{02}(T_2-T_0)-{\lambda }^E\min ({\rho }_0{\ell }_0\frac{dU}{dt},0).\\ C_1\frac{d{T}_1}{dt}={\tau }_0{\alpha }_{1}W+\sigma ({\epsilon }_0^d{T}_0^4-{\epsilon }_1{T}_1^4)-h_{01}(T_1-T_0)+{\eta }_1-{\lambda }^E\max ({\rho }_0{\ell }_0\frac{dU}{dt},0)\\ C_2\frac{d{T}_2}{dt}={\tau }_0{\alpha }_{2}W+\sigma ({\epsilon }_0^d{T}_0^4-{\epsilon }_2{T}_2^4)-h_{02}(T_2-T_0)+{\eta }_2-{\lambda }^E\max ({\rho }_0{\ell }_0\frac{dU}{dt},0)\\ \frac{dU}{dt}=(p{\mu }_1+q{\mu }_2)U_s\circ T_0-(p{\mu }_1+q{\mu }_2+\nu )U\\ \frac{d\vartheta }{dt}=\frac{2\pi }{Y{\sqrt{1-e^2}}^3}(1+e\cos \vartheta {)}^2.\end{array}$$

$${\alpha }_0^T=0.8,h_{02}=25,{\mu }_1=5.7\times 10^{-5},\nu =1.7\times 10^{-5}.$$

$$C_1=3.2\times 10^6,{\ell }_2=40,{\alpha }_0^T=0.88,h_{02}=30,{\mu }_1=2.9\times 10^{-6},\nu =10^{-5}.$$

$${\alpha }_0^T=0.85,h_{02}=35,{\mu }_1=5.7\times 10^{-5},\nu =2\times 10^{-5}.$$

$${\alpha }_0^T=0.84,h_{02}=35,{\mu }_1=4.6\times 10^{-5},\nu =2\times 10^{-5}$$

$${\alpha }_0^T=0.75,h_{02}=6,{\mu }_1=10^{-4},\nu =5\times 10^{-5},$$

$$\begin{cases}\dot{Q}_{1}=-(a+c)Q_{1}+d\,Q_{2}+\\ \qquad\qquad+\delta_{1}\big{[}\alpha\sin(\omega\,t)\sin(\Omega\,t)+\beta\cos(\omega\,t)\cos(\Omega\,t)+\gamma\sin(\omega\,t)\big{]}\\[5.0pt] \dot{Q}_{2}=c\,Q_{1}-(b+d)Q_{2}+\\ \qquad\qquad+\delta_{2}\big{[}\alpha\sin(\omega\,t)\sin(\Omega\,t)+\beta\cos(\omega\,t)\cos(\Omega\,t)+\gamma\sin(\omega\,t)\big{]}.\end{cases}$$

$$\begin{cases}\dot{Q}_{1}=cQ_{2}-(a+c)Q_{1}+\delta_{1}\Big{[}(\alpha+\beta)\cos(\Omega_{-}t)+(\beta-\alpha)\cos(\Omega_{+}t)+\gamma\sin(\omega\,t)\Big{]}\\[5.0pt] \dot{Q}_{2}=cQ_{1}-(b+c)Q_{2}+\delta_{2}\Big{[}(\alpha+\beta)\cos(\Omega_{-}t)+(\beta-\alpha)\cos(\Omega_{+}t)+\gamma\sin(\omega\,t)\Big{]},\end{cases}$$

Full text

The model takes into account the doubly-periodic forcing of the solar radiation (annual and daily), the laws of irradiance, the fact that the Earth has land and oceans with different thermic inertia, and the humidity of the air due to evaporation. This enables us to reproduce remarkable features of Earth’s climate such as lag of seasons, lag of noons, and asymmetric evolution of daily temperatures.

The model can easily be adapted to planets with non-terrestrial astronomic parameters. We conclude this article with an investigation of an Earth with eccentricity higher than real.

Key words. nonlinear dynamical systems, climate modelling, local climates, lag of seasons, earth and planetary climate

AMS subject classifications. 34C60, 37M05, 37N05

1 Introduction

Since the origin of modern meteorology, in the late 1800, researchers have based their models on a balance between incoming and outgoing radiative energy. Meteorology is in fact driven by such energy, and its main macroscopic indicator, the temperature, depends on the thermal inertia of the materials being irradiated and on the heat exchange among the different materials composing the surface of the planet. Relying on these ingredients, early meteorologists considered weather forecasting impossible [14].

Rapid local variations of temperature and other measurable quantities are due to motion of fluid masses and, excluding systematic effects due to Coriolis force and conformation of oceans bed and shores, such variations reasonably average out. It follows that the evolution of average macroscopic thermodynamic quantities is reasonably not as much influenced by such phenomena [15]. This is the main difference between meteorology and climate. For the reasons above, many aspects of climate can, and maybe should, be investigated disregarding meteorological models. In particular, many climatic effects can be exposed reconsidering Margules and Richardson’s basic models which, despite their simplicity, can give accurate climatic predictions.

In this article we focus on three main climatic effects. The first phenomenon is called lag of seasons. With this name one indicates the well known fact that the warmer days of the year take place some time after the days of maximal solar irradiance. The second phenomenon is the less celebrated phenomenon of lag of noons. With this name we indicate the fact that the warmer hours of the day take place some time after the hours of strongest solar irradiance. The third phenomenon is the asymmetries of temperatures. With this name we indicate the fact that daily temperatures rise much faster in the morning than they fall in the afternoon; this happens despite the fact that the forcing term is perfectly symmetric in shape, and an approach that takes into consideration only sun’s irradiation and Fourier’s law cannot reproduce this phenomenon. To obtain a realistic shape of daily temperatures one needs to introduce humidity in the picture.

The ultimate goal of this manuscript is to suggest a model that not only reproduces the three climatic effects above, but also predicts an extremely realistic evolution of average temperature and humidity. Only now accurate meteorological datasets of virtually any place in the world and for about 50 years are available. For this reason only now it is possible to compute a climatic unfolding of temperatures and other meteorological parameters (e.g. humidity) for any region on our planet.

To relate our work with the literature, we recall that in [3] the authors discuss a very simple mathematical model to explain the phenomenon of lag of seasons. Elementary mathematics proves that the long-term solution to the equation

[TABLE]

is

[TABLE]

where $\tau=\varphi/\omega$ and $\varphi=\arg(\lambda+i\omega)$. Equation (1) is an extremely basic model for the evolution of temperatures of a region $\mathcal{R}$ on the surface of a planet. In this model $T$ is the temperature of $\mathcal{R}$. The term $-\lambda T+\mu$ is a linearisation of the outgoing radiation from $\mathcal{R}$ while the forcing term $a\cos(\omega t)+b$ models the solar irradiance absorbed by $\mathcal{R}$. It clearly follows from (2) that the temperature $T$ has maxima and minima delayed with respect to the maxima and minima of solar irradiance, and the lag of these extreme temperatures is $\tau$.

The simple introduction of a forcing term containing two frequencies (daily rotation and annual revolution) is not rich enough to give realistic predictions of both lags (noons and seasons) [18]. To make such predictions one must increase the number of degrees of freedom. In fact, the first model whose solutions correctly predict both effects uses at least two different thermodynamic bodies with different thermic inertia, which correspond to a system with two degrees of freedom[6]. A quantitative analysis of the two lags is performed in Appendix 6. The asymmetry of temperatures is a further effect, and to be reproduced it requires the introduction of one more degree of freedom, that models the evolution of the absolute humidity of the air.

In the literature, models can be roughly divided into two categories: global circulation models (GCMs) [15] and energy balance models (EBMs) [3]. In GCMs land, oceans, and atmosphere are discretised into cells, and flows and energy transfer among cells are integrated over time; in EBMs the evolution of temperature is computed through low-dimensional systems, and the investigation is typically local or mediated along a parallel. GCMs can predict climate more accurately, but they require great effort to acquire data, to set up the simulation, and need large computing capacities. EBMs are possibly less accurate but require much less computational resources. EBMs have often been used to investigate climate under hypothetical variations of orbital and environmental parameters [26, 22]. Our work belongs to this second class of models. Our model could be used to investigate possible climate changes on Earth (e.g. greenhouse effect) and climate habitability of exoplanets in specific parts of their surface. In particular, at difference from classical one-dimensional EBMs [21], our approach is applicable when the revolution period and the rotation period are in 1:1 resonance (tidal locking) or other low-order resonance.

The outline of the work is the following. In Section 2 we give some geometric definitions and we write explicitly the expression of solar rays inclination. In Section 3 we recall the general laws of heat exchange and evaporation, and we write the evolution equations for temperature and humidity of a planet’s region. In Section 4 we numerically solve the equations for various regions on Earth, showing that our model well describes different types of climates (according to the Köppen climate classification [2]). In Section 5, we discuss the results and we indicate possible improvements and applications of the model.

2 The geometry of solar radiation

The motion and orientation in space of a region $\mathcal{R}$ on the surface of a planet is in good approximation due to the composition of the Keplerian revolution of the planet around its star and the rotation of the planet around its axis. The combination of such motions determines intensity and angle of the solar radiation responsible for the heating of the region. Disregarding all possible perturbations to this setting, the power of incoming solar radiation in $\mathcal{R}$ is hence completely determined by its exposition on the planet and the position of the planet in space.

2.1 Geometrical definitions

In our model the planet is assumed to be spheric. Its center of mass, following Kepler’s laws, revolves around the sun along an ellipse belonging to a plane called ecliptic plane. The planet also rotates uniformly around an invariable axis which makes a fixed angle $\gamma$, called obliquity, with respect to the normal of the ecliptic plane. The two points of the planet whose movement is not due to rotation are called North and South poles, and we agree that they are respectively at latitude $+90$ and $-90$ degrees (or $\pi/2$ and $-\pi/2$ radiants). The tropics are the two circles of points that have latitude $\pm\gamma$. We plan to describe the evolution of temperature in a certain region of the planet situated at a fixed latitude $\varphi$ and longitude $\psi$.

Astronomically speaking, significant instants are those in which the sun rays have local and global minimal (or maximal) distance from the zenit. Climatically speaking, significant instants are those in which the temperature has local and global maximum (or minimum). We hence give the definition of such events.

Definition 1

The solar solstices are the two instant in which the sun is at the zenit in one of the two tropics. The thermal solstices are the global extremes (maximum and minimum) of the temperature in a zone of the planet during the year.

Let us observe that at the solar solstice the following equivalent facts also take place:

• •

the projection of the terrestrial axis on the ecliptic plane is along the planet-sun line, and the pole in the same hemisphere of the zone is exposed to the sun;

• •

the sun is at the highest point when seen either from the North pole or from the South pole.

Definition 2

The solar noon is the instant in which the sun is at the local maximal height with respect to the horizon. The thermic noon is the moment in which the temperature is at a local maximum.

As anticipated in the Introduction, two remarkable phenomena take place on Earth: the lag between solar and thermal solstices and the the lag between solar and thermal noons.

Definition 3

The lag of seasons is the delay between the thermal and the solar solstice (summer and winter). The lag of noon is the delay between the thermal and the solar noon.

The astronomical special positions called aphelion and perihelion are unrelated to solstices. When the revolution of the planet around the star is not circular, we will have to keep into consideration the shift between summer/winter solstices and such aphelion/perihelion.

2.2 Inclination of solar rays

Let us consider a planet $P$ rotating around its sun $S$, and let $e_{1},e_{2},e_{3}$ be an orthonormal reference frame fixed with respect to the stars. The vector $e_{1}$ is parallel to the major semiaxes of the keplerian orbit of $P$ and is directed from $S$ to $P$ when $P$ is at the perihelion; the vector $e_{3}$ is normal to the ecliptic plane and is such that the rotation of $P$ around the sun is counterclockwise; the vector $e_{2}=e_{3}\times e_{1}$ completes the frame and is parallel to the minor semiaxes.

Following the classical description of keplerian motions, and supposing that at time $t_{0}=0$ the planet $P$ is located at the perihelion, the position of $P$ with respect to the sun is given in polar coordinates by the formulas

[TABLE]

[TABLE]

of reference frame $f_{1},f_{2},f_{3}$ attached to the rotating body with $f_{3}$ parallel to the axis of rotation is

[TABLE]

The versor connecting the planet to the sun is

[TABLE]

Since the region $\mathcal{R}$ at latitude $\varphi$ and longitude $\psi$ has normal to the surface

[TABLE]

it follows that

[TABLE]

This scalar product will be used in the following section, when writing the solar irradiance. In Table 1 the values of all parameters used in this discussion are indicated when the planet is the Earth.

3 The physics of heat transfer

The temperature of a region $\mathcal{R}$ on a planet is the result of a balance between the incoming radiation from the sun, the outgoing radiating energy, latent heats and heat exchanges within the system. We model the heat dynamics of such limited region $\mathcal{R}$ located at a certain latitude $\varphi$ and a certain longitude $\psi$. We disregard spatial diffusion and hence we use ordinary differential equations in which time is the independent variable. This is not a reasonable assumption when dealing with meteorology, but climatology deals with average evolution of the temperatures and the influence of neighbouring regions should average out. It follows that, to model experimental mean data, it is reasonable to suppose that the region $\mathcal{R}$ is physically isolated from the rest of the planet.

We restrict our study to the lowest part of the atmosphere and to the superficial layer of the planet’s surface. As we said in the Introduction, in order to reproduce lags, daily patterns, and more generally local climates, we consider three different homogeneous thermodynamic bodies, that in the case of Earth are air (which temperature we measure), land, and sea. To keep the model simple, we consider a unique mixed layer of air [19] for the lower atmosphere. We refer to the air layer using the index 0, to the land using the index 1, and to the ocean using the index 2.

This model represents the energy balance of the region $\mathcal{R}$ of the Earth’s surface which extension is reasonably of the order of $100$\mathrm{k}\mathrm{m}^{2}$. The real value of this surface plays no role in our investigation. In fact all quantities will be expressed “per unit surface”, and the units will always be divided per$\mathrm{m}$^{-2}$.

3.1 Solar irradiance

Approximating the sun to a black body, the solar irradiance flowing through a unit area perpendicular to the rays at distance $\rho$ from the sun is given by the Stefan-Boltzmann law

[TABLE]

Here $\sigma$ is the Stephan-Boltzmann constant ($\mathrm{J}\mathrm{s}^{-1}\mathrm{K}^{-4}\mathrm{m}^{-2}$), $R_{s}$ is the radius of the sun, $T_{s}$ is the temperature of the sun (see Table 2). In order to have the effective power received by a unit region $\mathcal{R}$ on the planet, we must multiply (5) by the scalar product (4). Considering the fact that during the night the contribution of the solar radiation is zero, the solar irradiance on $\mathcal{R}$ is

[TABLE]

Observe that this quantity is expressed in $\mathrm{J}\mathrm{s}^{-1}\mathrm{m}^{-2}$ and is the power of solar irradiance per unit area. When a light ray hits a body, its energy can be absorbed, transmitted or reflected. These three phenomena can be modelled introducing three parameters: absorbance $\alpha$, transmittance $\tau$, and reflectance $r$ such that $\alpha+\tau+r=1$. We mention here that in the literature the fraction of reflected radiation is commonly called albedo.

The solar rays cross the whole atmosphere, which absorbs a part of them. When the rays reach the surface, a part of them is absorbed by the superficial layer, another part is transmitted to a deeper underlying layer and a last part is reflected back to the atmosphere. Again, a part of this reflected radiation is absorbed, reflected or transmitted by the atmosphere. The layer of land and sea absorb all the incoming radiation but in a very different way; for this reason we must keep in mind that each region $\mathcal{R}$ is partly land and partly water. For this reason we introduce the main climatic parameter: a number $p\in[0,1]$ which represents the fraction of land and is referred to as solid fraction parameter. Its complementary parameter $q=1-p$ is the fraction of ocean.

The quantity of solar radiation absorbed by the three layers follows the laws

[TABLE]

The superscript $SR$ indicates that the contribution comes from Solar Radiation. The quantities $Q_{i}$ are expressed in $\mathrm{J}\,\mathrm{m}^{-2}$ and represent the heat quantity of the three thermodynamic bodies per unit area. The true amount of energy stored in such bodies can be obtained multiplying by the surface taken into consideration.

The parameters $\tau_{i},\alpha_{i},r_{i}=1-\alpha_{i}-\tau_{i},p,q=1-p$ are considered constants. We are aware that they actually are slightly variable, depending on the zenith distance of the sun, the atmosphere composition, the superficial temperature, and other factors. We will use their average value in the numeric integration. In our simulations, we have chosen $r_{1}=0.2$ for the reflectance of the land, which is a good approximation for Earth continents [29]. For other types of surface we can consider values of $r_{1}$ in the range $[0.1,\,0.4]$ [27]. The lowest values are appropriate for basaltic rocks or conifer forests, Sahara’s desert has $r_{1}\simeq 0.4$ [16], while grasslands have $r_{1}\simeq 0.2$ [19]. With respect to the ice, it has been documented a difference between ices over lands and over oceans [19, 12]. Therefore, following [27], we adopt $r_{1}=0.85$ and $r_{2}=0.62$ for ices over lands and ices over oceans respectively. We also suppose that all of the solar radiation not reflected by the surface is absorbed, giving $\alpha_{1}=1-r_{1}$, and $\alpha_{2}=1-r_{2}$. For the atmosphere the absorbance of solar radiation $\alpha_{0}$ is slightly variable [25], we assign to it the average values $0.25$. The transmittances $\tau_{i}$ are given by the relation $\tau_{i}=1-\alpha_{i}-r_{i}$. In Table 3 the values of all relevant parameters are listed.

3.2 Thermal radiation

All hot objects radiate with a Stefan-Boltzmann law. Unlike the sun, warm objects cannot be assumed to be black bodies and hence the power of emitted energy is $\varepsilon\sigma T^{4}$, where $\varepsilon$ is the emissivity of the body, a number in $[0,1]$ which depends on chemical and physical properties of the hot body. In this model the atmosphere will be assumed to radiate in two directions, down towards the earth with emissivity $\varepsilon_{0}^{d}$, and up towards outer space with emissivity $\varepsilon_{0}^{u}$. We also assume that $\varepsilon_{0}^{d}>\varepsilon_{0}^{u}$ because of lower density and temperature of the upper part of the atmosphere, and that all downward infrared radiation is absorbed by soil and water. We choose $\varepsilon_{0}^{d}=0.8$ for the radiation to the earth surface and $\varepsilon_{0}^{u}=0.45$ for the radiation to outer space.

The correct energy balance at our temperatures must include a parameter $\alpha_{0}^{T}$ to model the absorbance by the atmosphere of the radiation, called thermal radiation, emitted from Earth [25]. Unlike solar radiation, the spectrum of thermal radiation is mainly infrared, and $\alpha_{0}^{T}$ is much higher than $\alpha_{0}$. The value assigned to $\alpha_{0}^{T}$ is connected to the modelling of the greenhouse effect and it belongs to the interval $[0.8,0.95]$.

Summarizing, the power of energy transferred through thermal radiation between the thermodynamic bodies in $\mathcal{R}$ is

[TABLE]

The superscript TR stands for Thermal Radiation. For the thermal radiation, we consider these values of emissivity $\varepsilon_{1}^{soil}=0.94$ for soil, $\varepsilon_{1}^{sand}=0.75$ for deserts, $\varepsilon_{2}=0.96$ for oceans, and $\varepsilon_{1}^{ice}=\varepsilon_{2}^{ice}=0.85$ for ices over land and over ocean [1]. We suppose that the atmosphere absorbs most of the radiation emitted by the surface. All values are summarised in Table 4.

3.3 Conduction and convection

According to Fourier’s law, the rate at which two warm bodies exchange heat is proportional to the negative gradient of the temperature and to the area through which the heat flows. A similar law exists for convection, and is called Newton’s law of cooling. Altogether, if $T_{1}$ and $T_{2}$ are the temperatures of the two thermodynamic bodies, the heat flow $Q$ due to conduction and convection between them follows the law

[TABLE]

where $h$ is the cumulative heat transfer coefficient. In our model, the contributions of heat exchange due to conduction and convection are

[TABLE]

where $h_{ij}$ is the heat transfer coefficient among the two components labelled $i$ and $j$. In Table 5 the range for such coefficients are reported.

3.4 Geothermal heat.

In our model we take into consideration geothermal energy, that is heat coming from the mantle. There is a well defined region separating the mantle from the planet’s crust, called *Mohorovičić discontinuity *or Moho. Since the temperature of the mantle is much higher than the temperatures on the surface, we can assume that the geothermal heat flow is constant and we write

[TABLE]

The parameter $\eta_{i}$ is the power of energy conducted from the mantle to the body $i$ per unit area. Using experimental data from 20201 sites covering 62% of the Earth’s surface, Pollack et al. in [20] have obtained the values shown in Table 6.

The contribution of geothermal heat is between two and three order of magnitude lower than the contribution given by solar radiation. Its effect is hence feeble on Earth.

3.5 Evaporation.

The evaporation is a phenomenon that effects the absolute humidity of the air, and it depends on the wind speed and on the difference between the absolute humidity (the amount of kilograms of water vapour that a kilogram of dry air contains) and the saturation humidity of the air (the amount of kilograms of water vapour that a kilogram of dry air can contain at saturation). Saturated air cannot absorb water vapour, dry air does absorb vapour faster. We will assume that evaporation from land and sea is given by the law $\mu_{i}(U_{s}(T_{0})-U(t))$, with $\mu_{i}$ the rate of evaporation from the land and the sea expressed as a frequency ($\mathrm{s}^{-1}$). We assume the parameter $\mu_{2}$ to be variable depending on the wind speed. We consider the following approximation [8, 17]:

[TABLE]

The parameter $\mu_{1}$ is very variable, from very low values in the desert to very high values in tropical forests.

Measurements on Earth’s surface indicate that averaged wind speed is very different from one place to another, varying in $[0.17,\,20]\,$\mathrm{m}\mathrm{/}\mathrm{s} [[10](#bib.bib10)] at 10 $\mathrm{m}$ above the surface. When not directly accessible, we will use a reference value $\overline{v}_{i}=4\,$\mathrm{m}\mathrm{/}\mathrm{s}, which gives $\mu_{2}=5.7\times 10^{-5}\,$\mathrm{s}^{-1}$$.

Another factor that subtracts water vapour from the atmosphere is rain. The physical process that causes rain is condensation when the moist air rises to higher and colder strata of the atmosphere. We model this effect assuming a rate of rainfall proportional to the absolute humidity of the low atmosphere, and we call the coefficient of proportionality $\nu$, espressed in $\mathrm{s}^{-1}$. The values of $\nu$ can be computed knowing average rain precipitation in a year $\Pi$ (in meter of rain per square-meter), average humidity of the air $\overline{U}$ in kilograms of water vapour per kilogram of dry air. The parameter $\nu$ can be obtained using the formula

[TABLE]

We indicate with $\rho_{0}$ the density of air, with $\rho_{2}$ the density of water, and with $\ell_{0}$ the depth of the atmospheric layer. It turns out that reasonable values for $\nu$ are of the order of $10^{-5}$.

Summarizing, the law that regulates the evolution of humidity in the air is

[TABLE]

where

[TABLE]

is the humidity of saturation, a function of the air temperature $T$ whose values are the maximal amount of $\mathrm{kg}$ of water vapour that a $\mathrm{kg}$ of dry air can contain. This function has been obtained fitting well known values, its graph is represented in Figure 1.

3.6 Latent heat of evaporation and condensation.

As we have seen in last section, humidity plays a crucial role in the thermodynamics of the system under investigation. In fact, given a certain temperature and atmospheric composition, evaporation and condensation of water take place, depending on the difference between the absolute humidity and the saturation humidity. As we know, for each phase transition there is a latent heat, that is heat used for phase transition. During the process of evaporation part of the solar energy is used to change from liquid to vapour phase. That energy is not used to increase the temperature of the thermodynamic body. Therefore, if the mass of water undergoing evaporation per unit time and area is given by $\rho_{0}\ell_{0}dU/dt>0$, the related latent heat of evaporation is given by the formula

[TABLE]

where $\lambda$ is the specific latent heat for evaporation of water [23]. If $dU/dt<0$, the opposite process, called condensation, takes place. During this process, heat is released to the environment, with the same law as that for evaporation. In our system, the latent heat of condensation is released directly to the atmosphere, with the law

[TABLE]

Averaging during the year, it is known that heat exchanged through these processes amounts at about 25% of the solar irradiance [28]. To compare the magnitude of this process, heat transfer through convection amounts to about $5\%$ of the solar irradiance, and the energy absorbed directly by the atmosphere is between $18\%$ and $25\%$ of solar irradiance.

3.7 Thermal inertia.

Under the effect of heat transfers, the rate at which the temperature of a thermodynamic body change depends on its thermal capacity. In our case

[TABLE]

The parameters $C_{i}$ are the thermal capacities per unit surface ($\mathrm{J}\mathrm{K}^{-1}\mathrm{m}^{-2}$). For a body $i$ with density $\rho_{i}$ ($\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$), specific heat capacity $c_{i}$ ($\mathrm{J}\text{\,}{\mathrm{K}}^{-1}\text{\,}{\mathrm{kg}}^{-1}$), and depth $\ell_{i}$, the thermal capacity per unit surface is $C_{i}=\rho_{i}\,c_{i}\,\ell_{i}$. We will assume the thermal capacities constant for all thermodynamic bodies except for the air. This is justified by the fact that the thermal capacity of the air depends on its content of water vapour. Recalling that $U$ is the absolute humidity of the air, measured in $\mathrm{kg}$ of water per $\mathrm{kg}$ of air, we will assume that the heat capacity of the air is

[TABLE]

where $c_{0}^{d}$ is the specific heat capacity of dry air, and $c_{0}^{v}$ is the specific heat capacity of water vapour, $\rho_{0}^{d}$ is the density of dry air, and $\ell_{0}$ is the effective depth of the layer of air. The new independent variable $U$ here introduced will in turn depend, via a differential equation, from the temperature of the air. The specific heat of dry air is $711.28\,$\mathrm{J}\text{,}{\mathrm{kg}}^{-1}\text{,}{\mathrm{K}}^{-1}, the specific heat of water vapour is $2050.16\,$\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}\text{\,}{\mathrm{K}}^{-1}. In our model we consider a layer of lower atmosphere $\ell_{0}=400\,$\mathrm{m}$$.

Following [19, 29], we choose $\ell_{1}\in[0.3,\,0.5]\,$\mathrm{m} for soil, $\ell_{2}\in[40,\,60]\,$\mathrm{m} for oceans. The thermal characteristics of land and water differ from region to region. For this reason in different cases we use different heat capacities. For details on such values see [27, 11, 4, 7]. Using the arguments above one obtains the values shown in Table 8.

3.8 Final system.

Summarizing equations (3), (4), (6), (7), (8), (10), (11), (17), (13), (15), (16), with some minimal algebra, the dynamical system that models the temperature evolution of $\mathcal{R}$ is modelled by the evolution of the 4 independent variables $T_{0},T_{1},T_{2},U$, and the variable $\vartheta$, whose evolution is fixed for planet Earth,

[TABLE]

4 Numerical analysis

In this section we justify the final choice of the parameters depending on the choice of region $\mathcal{R}$ on Earth, and we run the simulation for various types of climates. We then compare the numerical results with the real average temperatures. The numerics, the acquisition of real temperatures, and their manipulation have been done using the software Mathematica Wolfram Research Inc. In particular WeatherData[] allowed us to acquire the dataseries of temperatures and humidity with respect to coordinated universal time (UTC) from a variety of weather stations in regions with different climates. We averaged the temperatures at any given time of the year over a period of 47 years (from 1973 to 2019).

We consider 5 regions: Hilo–Hawaii, Kufra–Lybia, Catania–Italy, Lincoln–USA,Vostok–Antarctica. Each region belongs to one of the Köppen climate zones [2]: Tropical (A), Arid (B), Temperate (C), Continental (D), and Polar (E). In the following subsections we choose the local parameters and we superimpose mean temperature and humidity in the chosen region with temperature and humidity obtained with our model. In the yearly plots we indicate with dashed lines the solstices and equinoxes; in the daily plots we indicate midday and midnights.

4.1 Tropical climate: Hilo, Hawaii

Hilo belongs to a region with Tropical, Rainforest Köppen climate (AF type). It is situated at latitude $19.72$ and longitude $-155.05$. Belonging to an Hawaiian island, we choose $p=0.1$. The presence of forest makes it reasonable to choose an higher value for the land’s thermal capacity $C_{1}=1.7\times 10^{6}$\mathrm{J}\mathrm{K}^{-1}\mathrm{m}^{-2}, while considering $\ell_{2}=50$\mathrm{m} for oceans gives $C_{2}=2.1\times 10^{8}$\mathrm{J}\mathrm{K}^{-1}\mathrm{m}^{-2}$$. We also consider the following values for other location-dependent parameters:

[TABLE]

In Figure 2 are represented the computed evolution of temperature (red) and humidity (blue) of the air and the real averaged temperatures and humidities (black) from 1973 to 2019.

4.2 Arid climate: Kufra, Libya

Kufra belongs to the eastern part of Sahara with Arid, Hot Desert (BWH type) Köppen climate. It is situated at latitude 24.18 and longitude 23.31. Being in a desert, water has almost no influence and we hence have chosen $p=0.9$. We recall the choices

[TABLE]

All other parameters are in the Tables, and we have used the coefficients for sand. In Figure 3 are represented the computed evolution of temperature (red) and humidity (blue) of the air and the real averaged temperatures and humidities (black) from 1973 to 2019.

4.3 Temperate climate: Catania, Italy

Catania is one of the cities on the Mediterranean Sea with Temperate, Hot-summer, Mediterranean Köppen climate (CSA type). It is situated at latitude 37.47 and longitude $15.05$. Given its location, we choose $p=0.6$. Considering that the top layer of Mediterranean sea mix to a depth of up to $40$\mathrm{m}$=\ell_{2}$, we consider

[TABLE]

In Figure 4 are represented the computed evolution of temperature (red) and humidity (blue) of the air and the real averaged temperatures and humidities (black) from 1973 to 2019.

4.4 Continental climate: Lincoln, USA

Lincoln belongs to the central USA, a region with Continental, Hot-summer, Humid Köppen climate (DFA type). It is situated at latitude $40.85$ and longitude $-96.75$. For its location, we choose $p=0.8$, $\ell_{2}=40$\mathrm{m}$$. We adopt

[TABLE]

In Figure 5 are represented the computed evolution of temperature (red) and humidity (blue) of the air and the real averaged temperatures and humidities (black) from 1973 to 2019.

4.5 Polar climate: Vostok, Antarctica

Vostok is a weather station close to a lake in Antartica, it is located almost at the South Pole and it has Polar, Ice cap Köppen climate (EF type). This region is situated at latitude $-78.45$ and longitude $106.87$ and is always covered with ice and snow, living in eternal winter. The thermal inertia of the ice cap is very high and so, even if located on land, we have chosen $p=0.45$,

[TABLE]

All other parameters are in the Tables, and we have used the coefficients for ice over land and over water.

In Figure 6 are represented the computed evolution of temperature (red) and humidity (blue) of the air and the real averaged temperatures and humidities (black) from 1973 to 2019.

4.6 More eccentric cases

It is well known that the for solar system the major semiaxis of the Earth’s orbit is stable under perturbations while the stability of the full set of orbital parameters is still much discussed in modern times [13]. In our model the orbital parameters can be easily changed to model the temperatures in an Earth-like planet. The small eccentricity of the orbits in the solar system are well known to be non-generic [9]. In the following plots we investigate the temperatures that Catania would have if the eccentricity of Earth was $e=0.2$ or $e=0.5$, and we compare the same effect on Sydney, a city in the southern hemisphere. We recall that, because of Earth’s orientation of the rotation axis, during the summer of the northern hemisphere the Earth is at the aphelion, while during the summer of the southern hemisphere the Earth is at the perihelion. It follows that the effect of a change in eccentricity is mild in Catania (see Figure 7 top) and severe in Sydney (see Figure 7 bottom). Let us note however that the precession of the equinoxes would switch the situation every half Platonic year (12886 years).

These speculations are particularly interesting for their applications on exoplanets, where the suitability of temperatures to host life is a fundamental issue [5, 24].

5 Conclusions

In this article, starting from basic physical laws, observations on the geometry of solar systems, and knowledge of the structure of a given planet, we design a model for local climate. The parameters involved in this dynamical model are mostly given by experimental experience. Our investigation is restricted to Earth, for which we can reproduce climatic phenomena like lag of seasons, lag of noons, asymmetry of daily temperatures, evolution of temperatures related a variety of Köppen climatic zones.

The temperatures computed solving the equations are reasonably similar to the real ones. In particular, as the real one, display lags and asymmetric evolution. They could be better fitted with a more detailed choice of the parameters or adding some other phenomenon to the equations. The annual excursion is very reasonable in all models, the daily temperature excursions tend to be slightly smaller than real (see Figures 3–5). Despite the fact that humidity creates an asymmetric raise and fall of temperatures, the raise of real temperatures in the morning is faster than simulated ones. Simulated humidities tend to be constant during the year, while in some regions real humidities have smaller values in the summer than in the winter; the daily excursions are much more accurate.

The model still displays some criticalities and it can be improved in many ways. In particular we indicate the following issues:

the model requires the inclusion of some water also when dealing with desert or ice-caps, because the stabilising effect of water is necessary to avoid high annual temperature excursions; 2. 2.

we only consider the lowest part of the atmosphere. 3. 3.

spatial diffusion has been disregarded; 4. 4.

the equations that model humidity is not completely satisfactory, it probably should take into account other factors;

The first and second issues could be dealt with by adding other layers, one below the soil and one above the lower atmosphere. This would grant a correct annual excursion without compromising the daily one.

The spatial diffusion has been intentionally excluded to keep the model as simple as possible. The introduction of diffusion completely changes the approach, forcing a discretisation of the surface of the planet and the creation of a GCM which requires a detailed description of the planet surface and a large computational effort.

The evaporation depends on wind velocity, and probably non-constant wind speed should be taken into account, as well as seasonal rainfall ratio variation. We have not made a deep investigation on this facts, and we do not propose solutions.

The investigation is suitable to applications to exoplanets. Some astronomical parameters of exoplanets are known, but the choice of most other parameters is a delicate issue and will be subject for future works. In particular we think that the model could be most useful in the investigation of habitable and tidally locked planets.

Acknowledgements

AG wishes to thank Giancarlo Benettin for his suggestion of considering two thermodynamic bodies in the model. GDB thanks Alberto Chiavetta for his help in understanding the physics of the system. Both authors thank Paolo Falsaperla for enlightening discussions on energy balance and an anonymous referee for his critical reading of a first version of the article.

AG is supported by the group GNFM of INdAM, by the grant No. PTRDMI-53722122113 of the University of Catania, and by Grant 2017YBKNCE of the national project PRIN of MIUR. GDB acknowledges support from Scuola Superiore di Catania.

6 Appendix: the mathematical essence

In this appendix we make a cumbersome mathematical analysis of the fundamental reason that justifies the double-lag phenomenon. To model the temperature evolution of two thermodynamic bodies driven by a doubly-periodic forcing term, we consider a system of two differential equations:

[TABLE]

We have discussed in Section 3 how this system models the temperature evolution of two different thermodynamic bodies in a zone of a planet. The only difference with equations (3.8) lays on the fact that the longitude is absent and the exchange of heat is not mediated by a layer of air. The two bodies are irradiated by solar rays modulated by two frequencies $\omega$ and $\Omega$ that are respectively $2\pi$ times the reciprocal of a year and $2\pi$ times the reciprocal of a day. The terms $a\,Q_{1}$ and $b\,Q_{2}$ model the heat flow from the bodies to space, the terms $c\,Q_{1}$ and $d\,Q_{2}$ model the rate of heat exchange among the two bodies.

Using Prostaferesi-Werner formulaes one can rewrite the equations as

[TABLE]

with $\Omega_{-}=\Omega-\omega$ and $\Omega_{+}=\Omega+\omega$. The homogeneous linear system associated to the equations has a stable node for every choice of parameters. In fact the determinant of the associated matrix is

[TABLE]

and the discriminant $tr^{2}-4\,det$ is

[TABLE]

With a linear change of variables of matrix $(S_{1},S_{2})=P(Q_{1},Q_{2})$ the system becomes

[TABLE]

where the vector $(\chi_{1},\chi_{2})=P(\delta_{1},\delta_{2})$, $P$ is the matrix of change of basis, $S=PQ$, and $\lambda_{1},\lambda_{2}$ are the two eigenvalues of the linear system. The actual expression of the coefficients $\lambda_{1}$, $\lambda_{2}$, $\chi_{1}$, $\chi_{2}$ is irrelevant for our purposes. What is important is that the asymptotic solutions to these equations have the form

[TABLE]

with

[TABLE]

for $i=1,2$. Turning back to the temperatures $T_{1},T_{2}$ one has

[TABLE]

The term

[TABLE]

is responsible of the yearly delay, that can be estimated with the following algebraic steps:

[TABLE]

The lag of seasons for the bodies 1 and 2 are the two components of the vector

[TABLE]

The lag of noon is more delicate. In fact the delay can be estimated only if the ratio $\Omega/\omega$ is large (as it happens on Earth). In such case the evolution of temperatures is the sum of two terms:

[TABLE]

If the ratio $\Omega/\omega$ is large, then

[TABLE]

It follows that the lag of noon of the two bodies is given by

[TABLE]

More precisely, one has that the solutions to the equation (18) are

[TABLE]