Milestone 2 - Solar forcing

  1. Important definitions
  2. Effect of Earth's distance to the sun
  3. Effect of Earth's axial tilt
    1. Latitudes where there is no sunrise (winter)
    2. Latitudes where there is no sunset (summer)
    3. Latitudes where there is daily sunrise and sunset
  4. Bonus: Computation of the true longitude λ\lambda

So far, we have made the assumption that the solar forcing term remains constant across both space and time. However, this assumption is again too simple, as it overlooks the fact that the Earth's position and orientation in relation to the sun changes throughout the year, leading to a variation in the amount of incoming insolation.

By incorporating these effects, we can significantly improve the modeling of our solar energy source term. The two primary factors that contribute to the variability of solar radiation received by the Earth are as follows:

  1. Distance from the sun: Rather than a perfect circle, the Earth's orbit around the sun is an ellipse. Consequently, the distance between the Earth and the sun fluctuates throughout the year, leading to a change in the total amount of solar radiation that reaches the Earth.

  2. Tilted axis of rotation: Due to the Earth's tilted axis of rotation, different parts of the planet receive varying amounts of solar radiation as it orbits the sun (depending on the latitude). This phenomenon gives rise to the changing seasons.

The derivations presented in this section can be found in the following references:

Berger, A. (1978). Long-term variations of daily insolation and Quaternary climatic changes. Journal of Atmospheric Sciences, 35(12), 2362-2367.

Berger, A., Loutre, M. F., & Tricot, C. (1993). Insolation and Earth's orbital periods. Journal of Geophysical Research: Atmospheres, 98(D6), 10341-10362.

Important definitions

Before we delve into calculating the impact of the Earth-sun distance and Earth's axis tilting on the solar radiation on Earth, let us define some relevant concepts (see figure above for reference):

  • Perihelion: The point in Earth's orbit where it is closest to the Sun. Currently, the perihelion occurs around January 4th.

  • Aphelion: The point in Earth's orbit where it is farthest from the Sun. Currently, the aphelion occurs around July 6th.

  • Ecliptic: The plane of Earth's orbit around the Sun.

  • Northern summer solstice: The point in Earth's orbit where the longest day of the year occurs in the Northern Hemisphere and the shortest day of the year occurs in the Southern Hemisphere. At this point, the plane formed by the rotation axis of Earth and the line that connects Earth and the Sun is perpendicular to the ecliptic. Currently, the Northern summer solstice occurs around June 21st.

  • Northern winter solstice: The point in Earth's orbit where the shortest day of the year occurs in the Northern Hemisphere and the longest day of the year occurs in the Southern Hemisphere. At this point, the plane formed by the rotation axis of Earth and the line that connects Earth and the Sun is perpendicular to the ecliptic. Currently, the Northern winter solstice occurs around December 21st.

  • Vernal equinox: The point in Earth's orbit where the line that connects Earth and the Sun aligns with the equatorial plane of Earth during the transition from winter to summer in the Northern Hemisphere. This is also the moment when Earth's rotation axis is directly perpendicular to the Sun-Earth line.

  • Northern autumnal equinox: The point in Earth's orbit where the line that connects Earth and the Sun aligns with the equatorial plane of Earth during the transition from summer to winter in the Northern Hemisphere. This is also the moment when Earth's rotation axis is directly perpendicular to the Sun-Earth line.

  • True longitude of Earth (λ\lambda): Position of Earth at a given time as measured from the vernal equinox. The position λ=0\lambda = 0 is considered the beginning of an astronomical year.

To compute the distance from the sun and Earth's axis tilt angle, we need to consider three parameters that vary over time as a result of the gravitational interactions between Earth and other celestial bodies within the solar system. We call these parameters the orbital parameters and define them as

  1. Eccentricity (ee): This parameter determines the degree to which Earth's orbit around the sun deviates from a perfect circle. It is computed as

e=1−b2a2, e = \sqrt{1-\frac{b^2}{a^2}},

where aa represents the semi-major axis, and bb represents the semi-minor axis of the ellipse (refer to the figure). Earth's eccentricity undergoes periodic changes over hundreds of thousands of years.

  1. Obliquity (ϵ\epsilon): This parameter, also referred to as axial tilt, represents the angle between Earth's rotational axis and its orbital axis. The obliquity angle is equivalent to the angle between the equatorial plane and the orbital plane. The obliquity angle varies periodically over approximately τϵ≈20,000\tau_{\epsilon} \approx 20,000 years.

  2. Precession distance (ω~\tilde \omega): Earth's axis rotates with a period of about 40,00040,000 years, causing the seasons to shift in position within the elliptical orbit of Earth. This phenomenon is known as precession or spin. We define the precession distance as the angle between the aphelion and the vernal equinox.

NOTE: In the paper of Berger, the precession distance (ω~\tilde \omega) is sometimes defined from the perihelion and sometimes from the aphelion. We use the aphelion in our implementation as it is the one obtained with equation (6) from Berger.

In the following sections, we will assume that the eccentricity, obliquity and precession distance are constant to derive the total solar irradiance. This is a reasonable assumption for climate simulations in which the orbital parameters are updated at the time interval Δtorb≪τϵ\Delta t_{\text{orb}} \ll \tau_{\epsilon}.

Effect of Earth's distance to the sun

The inverse-square law states that the total radiation that is received from a source is inversely proportional to the square of the distance from the source.

As a result, the total radiation that is received at the top layer of Earth's atmosphere is

Sr(t)=r02r2S0, S_r(t) = \frac{r_0^2}{r^2} S_0,

where rr is the (time-dependent) Earth-sun distance, r0≈ar_0 \approx a is roughly the mean Earth-sun distance (defined as one anstronomical unit), and S0S_0 is the so-called solar constant, the mean solar elecromagnetic radiation per unit area measured on a surface perpendicular to the solar rays at a distance r0r_0.

Since Earth's orbit is an ellipse and the sun sits in one focus, we can calulate the Earth-sun distance as

r=a(1−e2)1−ecos⁡(ν), r = \frac{a (1-e^2)}{1 - e \cos (\nu)},

where the position (angle) of Earth with respect to the aphelion is

ν=λ−ω~ \nu = \lambda - \tilde \omega

The radiation at the top of the atmosphere is then computed as

Sr(t)=(1−ecos⁡(λ−ω~))2(1−e2)2S0. S_r(t) %= \frac{(1-e \cos (\nu))^2}{(1 - e^2)^2} S_0 = \frac{(1-e \cos (\lambda - \tilde \omega))^2}{(1 - e^2)^2} S_0.

Effect of Earth's axial tilt

Due to Earth's axial tilt, the solar irradiance depends on the latitude and the time.

We call the angle between the Earth-sun line and the equatorial plane of Earth declination angle and note it with δ\delta. This angle changes throughout the year, as the Earth revolves around the sun (see figure):

In general, the declination angle can be computed as

sin⁡(δ)=sin⁡(λ)sin⁡(ϵ). \sin (\delta) = \sin (\lambda) \sin (\epsilon).

Depending on the declination angle and the latitude, we identify three regions of interest:

Latitudes where there is no sunrise (winter)

There is no incoming solar radiation for latitudes that fulfill

∣θ∣+∣δ∣≥π2 with θδ<0, |\theta| + |\delta| \ge \frac{\pi}{2} \,\, \text{with} \,\, \theta \delta < 0,

or equivalently

z0=−tan⁡(θ)tan⁡(δ)≥1. z_0 = -\tan(\theta)\tan(\delta) \ge 1.

For these latitudes, the insolation is simply

S(θ,t)=0. S(\theta,t) = 0.

Latitudes where there is no sunset (summer)

Some regions of Earth are exposed to incoming solar radiation throughout the entire day during part of the year. These regions fulfill

∣θ∣+∣δ∣≥π2 with θδ>0, |\theta| + |\delta| \ge \frac{\pi}{2} \,\, \text{with} \,\, \theta \delta > 0,

or equivalently

z0=−tan⁡(θ)tan⁡(δ)≤−1. z_0 = -\tan(\theta)\tan(\delta) \le -1.

The calculation of the effective insolation is more involved, as it requires the use of spherical astronomy. The result reads as

S(θ,t)=Sr(t)sin⁡(θ)sin⁡(δ). S(\theta,t) = S_r(t) \sin(\theta) \sin(\delta).

Latitudes where there is daily sunrise and sunset

There is daily sunrise and sunset at latitudes that fulfill the condition

−(π2−∣δ∣)<θ<π2−∣δ∣. - \left( \frac{\pi}{2} - |\delta| \right) < \theta < \frac{\pi}{2} - |\delta|.

For these latitudes, the calculation of the insolation also requires the use of spherical astronomy. The result reads as

S(θ,t)=Sr(t)1π(H0sin⁡(θ)sin⁡(δ)+cos⁡(θ)cos⁡(δ)sin⁡(H0)), S(\theta,t) = S_r(t) \frac{1}{\pi} (H_0 \sin(\theta) \sin(\delta)+\cos(\theta) \cos(\delta) \sin(H_0)),

with the absolute angle of the sun at sunrise and sunset

H0=arccos⁡(z0). H_0 = \arccos (z_0).

In summary, the insolation is given by the expression

S(θ,t)={0if ∣θ∣+∣δ∣≥π2 with θδ<0,S0ρ(t)sin⁡(θ)sin⁡(δ)if ∣θ∣+∣δ∣≥π2 with θδ>0,S0ρ(t)π(H0sin⁡(θ)sin⁡(δ)+cos⁡(θ)cos⁡(δ)sin⁡(H0))if −(π2−∣δ∣)<θ<π2−∣δ∣, S(\theta,t) = \begin{cases} 0 & \text{if}~|\theta| + |\delta| \ge \frac{\pi}{2} \,\, \text{with} \,\, \theta \delta < 0, \\ S_0 \rho(t)\sin(\theta) \sin(\delta) & \text{if}~|\theta| + |\delta| \ge \frac{\pi}{2} \,\, \text{with} \,\, \theta \delta > 0, \\ \frac{S_0 \rho(t)}{\pi} (H_0 \sin(\theta) \sin(\delta)+\cos(\theta) \cos(\delta) \sin(H_0)) & \text{if}~- \left( \frac{\pi}{2} - |\delta| \right) < \theta < \frac{\pi}{2} - |\delta|, \end{cases}

with the normalized distance between Earth and sun squared

ρ(t)=(1−ecos⁡(ν(t))1−e2)2. \rho(t) = \left(\frac{1 - e \cos (\nu(t))}{1-e^2}\right)^2.

The total solar forcing is given by

Ssol(θ,φ,t)=S(θ,t)(1−α(θ,φ)). S_{sol}\left(\theta,\varphi,t\right) = S\left(\theta,t\right) (1 - \alpha \left(\theta,\varphi\right)).

Bonus: Computation of the true longitude λ\lambda

Kepler's second law states that a line segment joining a planet and the Sun sweeps out equal areas during equal intervals of time. Therefore, Earth moves faster when it is closer to the sun. Is a short time interval dt\mathrm{d} t, the area swept is equal to the area of a triangle with base rdνr \mathrm{d} \nu and height rr. Therefore, we can write Kepler's second law as

dAdt=πabT=r22dνdt, \frac{\mathrm{d} A}{\mathrm{d} t} = \frac{\pi a b}{T} = \frac{r^2}{2} \frac{\mathrm{d} \nu}{\mathrm{d} t},

where πab\pi a b is the total area of a ellipse, which is swept in a complete period TT.

Assuming that the time is measured in years and replacing the definition of the Earth-sun distance (3) and the eccentricity (1) into (19) we obtain

dλdt=2π(1−e2)3/2(1−ecos⁡(λ−ω~))2. \frac{\mathrm{d} \lambda}{\mathrm{d} t} = \frac{2\pi}{(1 - e^2)^{3/2}} (1 - e \cos (\lambda - \tilde \omega))^2.
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.