Perpetual Winter

Winter is my favorite season. Refreshing rain, cold, comfortable clothes, clouds, and much less sweat, all make it superior to the hot summer. Why, oh why, can’t it be winter forever? It turns out, that if Earth’s astronomical characteristics were a little different, we could have had an eternal winter (at least in some places of the world).
It is more or less common knowledge, that the change of seasons and differences between the hemispheres is caused by Earth’s axial tilt, which is about 23.30 relative to the ecliptic (the plane of revolution around the sun).

Describes what the tilt of the Earth means

Generally, the further away from the equator a place on Earth is, the less energy it receives from the sun; this is because the ratio between its area and its cross section perpendicular to the ecliptic is larger than for lower latitudes. In the above picture, suppose that the sun is directly left of Earth. In this case, the southern hemisphere will occupy more of the area facing the sun, resulting in summer. Furthermore, parts of the northern half of the globe are continually hidden from the sun, creating a long night and a cold, harsh, winter.
The problem is, that once Earth gets to the other side of the sun, the exact same phenomena happen, but with the hemispheres reversed; once the sun is to the right of Earth, its tilt makes the northern hemisphere face more towards the sun, getting more sunlight and causing an awful summer. What we want, therefore, is for the tilt to change along with time. If, after half a year, Earth changes its tilt, so that when it gets to the other side of the sun, the northern hemisphere will still be on the farther and darker side, facing away from the sun, then we will have achieved perpetual winter. Of course, this also means perpetual summer for the southern hemisphere, but I suppose they will get used to it.

Shows what axial precession is

This phenomenon – the tilt changing its angle – is called axial precession. It is interesting to note that Earth does exhibit such movement, and its axis completes a single revolution once per about 25,800 years. Normally this is not noticeable, but be aware, that in several hundred years, Polaris will not be the Northern star any more. However, we want this effect to synchronize itself with Earth’s orbit around the sun. If this were so, then Earth’s orientation towards the sun would be “locked” – the same hemisphere would always point towards it.
This claim may seem rather intuitive, but it helps to visualise it. Take a small ball in your hand and draw a mark on the ball. Move only your arm in a circular path around a fixed point in space – this is Earth’s orbit around the sun. Now, rotate your hand as you move your arm, so that the mark on the ball always faces that fixed point in space – this is how we would like Earth’s tilt to act – axial precession.
We will now show mathematically that this effect generates the wanted result. This is a bit technical, but its worth it. Our main goal will be to see how much energy any given place on Earth receives when the axis undergoes precession. This will be done by looking at how much of the energy from the sun actually reaches a given location. We will use the following picture as our reference:

Shows the various rotations and movements of Earth, excluding nutation

With the x axis pointing at the sun, and the z axis pointing upwards, perpendicular to the ecliptic.
Let us consider what an observer who is orbiting the sun sees. He stands in the center of the world, at the very core of the planet, and circles about the sun just like Earth does; however, he does not rotate about himself at all, but his axes stay in a fixed orientation, the same as given in the picture. From his point of view, he is stationary, and the sun revolves around him. The vector pointing in the direction of the sun, as a function of time, is:

\mathbf{\hat{r}_{patch} }= \left( \begin{array}{c}\cos{\omega_{y}t} \\ \sin{\omega_{y}t} \\ 0 \end{array} \right)

where ωy is the angular velocity of Earth’s orbit around the sun (ωy = 2π / 1 year).
Let us take take a certain “tile”, a “patch” on Earth, as shown in the picture, and calculate its position vector. We will start with an Earth that is not moving and is untilted, and slowly build up its astronomical characteristics. Any position on Earth can be described by two angles: latitude, and longitude. For our purposes, as you will see, we may choose any longitude we wish without losing generality, and so choose 00. Hence, for a given latitude φ, the position vector of a patch is:

\mathbf{\hat{r}_{patch} }= \left( \begin{array}{c}0 \\ \cos{\varphi} \\ \sin{\varphi} \end{array} \right)

Now, let us consider Earth’s rotation around itself. We have assumed for now that it is untilted, so the axis of rotation and the z axis overlap. Thus, the rotation is simply around the z axis. Given our choice of axial orientation, as seen in the picture, a counter-clockwise rotation looks like:

\begin{array}{c} x_{new} = x_{old}\cos{\omega_{d}t} - y_{old}\sin{\omega_{d}t}  \\ y_{new} = y_{old}\cos{\omega_{d}t} + x_{old}\sin{\omega_{d}t} \\ z_{new} = z_{old} \end{array}

Where ωd is the angular velocity of Earth’s rotation around itself (ωd = 2π / 1 day). Applying this to our patch vector, we get

\mathbf{\hat{r}_{patch}} = \left( \begin{array}{c}-\cos{\varphi}\sin{\omega_{d}t} \\ \cos{\varphi} \cos{\omega_{d}t}\\ \sin{\varphi}\end{array}  \right)

This is the position of the patch, if Earth were untilted. Looking at our picture, we see that in order to tilt it, we have to rotate it an angle η around the observer’s y axis. This will give us the position of a path on Earth as it is today. Based on our axial orientation, the rotation is:

\begin{array}{c} x_{new} = x_{old}\cos{\eta} - z_{old}\sin{\eta}  \\ y_{new} = y_{old} \\ z_{new} = z_{old}\cos{\eta} + x_{old}\sin{\eta}  \end{array}

Applying this to our patch vector, we get

\mathbf{\hat{r}_{patch}} = \left( \begin{array}{c} -\cos{\varphi}\sin{\omega_{d}t} - \sin{\varphi}\sin{\eta} \\ \cos{\varphi} \cos{\omega_{d}t} \\ \sin{\varphi}\cos{\eta} - \cos{\varphi} \sin{\omega_{d}t}\sin{\eta}\end{array} \right)

Now all we have to do is add precession, and we know the final position of the patch as a function of time. Precession is once again, rotation around the z axis. You may be wondering, “wait, we rotated around the z axis when we considered Earth’s daily rotation around itself. Won’t these two just combine?” The answer is no; when rotating, the order in which you apply your operations matters. In this case, Earth’s axis and the observer’s z axis no longer align, because we already took into consideration Earth’s tilt. We will now apply the rotation, but around the observer’s z axis. The direction of rotation is the same as Earth’s rotation around the sun, so that later on we will be able to match it. We get (breathe deep):

\mathbf{\hat{r}_{patch}} = \left( \begin{array}{c} ( -\cos{\varphi}\sin{\omega_{d}t} - \sin{\varphi}\sin{\eta} ) \cos{\omega_{p}t} - \cos{\varphi}\cos{\omega_{d}t}\sin{\omega_{p}t} \\ \cos{\varphi} \cos{\omega_{d}t}\cos{\omega_{p}t} + ( -\cos{\varphi}\sin{\omega_{d}t} - \sin{\varphi}\sin{\eta} ) \sin{\omega_{p}t} \\ \sin{\varphi}\cos{\eta} - \cos{\varphi} \sin{\omega_{d}t}\sin{\eta}\end{array} \right)

Where ωp is the angular velocity of the precession.
We now have the precise position of a certain patch on Earth. As we have previously said, the amount of energy that a certain patch gets depends on the ratio between its surface area and its cross section. In other words, it depends on the angle between its normal and light coming in from the sun. Let j = S ⋅ rsun, where S is the solar constant – the energy density coming from the sun. Its current value is about 1368 [W/m2]. In this case, the total energy per square meter that a patch receives is dΦ = jrpatch. We note, however, that this dot product can generate negative results. This occurs when the cosine of the angle between the patch and the sun is negative. The physical interpretation, is that the patch is hidden, obscured from the sun – if the product is negative, then it is currently nighttime for that patch. In this case, it should receive no energy at all. We shall refine our equation then:

d\Phi = \left\{ \begin{array}{rcl} \mathbf{j}\cdot\mathbf{\hat{r}_{patch}}  & \mbox{if} & \mathbf{j}\cdot\mathbf{\hat{r}_{patch}} > 0 \\ 0 & \mbox{if} & \mathbf{j}\cdot\mathbf{\hat{r}_{patch}} \leq  0 \end{array}\right.

This gives the energy density for a given time, t. We want, however, to see the average density over a period of time. Then, we would be able to look at how much energy arrived every single day for a year, and compare the results for no precession, and for precession with a yearly period. This is also what lets us ignore the longitude, as mentioned earlier – if we take the average energy over time spans which are integer multiples of one day, there is no difference in longitude. Our final result to calculate, then, is:

\overline{\Phi} = \frac{1}{t_{2} - t_{1}}\int_{t_{1}}^{t_{2}} d\Phi dt

I do not know of any way to calculate this integral analytically. However, it is very easy to calculate it with a rather cheap numerical integration method, say the Trapezoidal Rule. For this purpose, a small python script was written, attached below (not intended to be robust or extensible).
So, how does the model fare? The graph given here depicts the energy density per day for a latitude of 31.50 (The Mediterranean’s latitude), with no precession. This reflects today’s status.

Energy density per day, no Earth precession, 31.5 latitude

The initial season, just like in the picture, is winter, hence the incident energy is very low. Half a year later, it is the peak of summer, and it is much hotter. The energy density almost reaches 500 [W/m2], which is quite impressive, as we might think the maximum is about 700 – half the time it is day, half the time it is night. Adding a precession period equal to exactly one year generates the following:

Energy density per day for yearly precession, 31.5 latitude

Just as expected! We get an eternal winter, with very little energy coming from the sun, all year long!

The script, although it is rather crude and simple, actually does a wonderful job of showing all sorts of nice effects. Reversing the latitude to -31.50 gives again a constant energy density, but at a value more than twice as large – meaning an eternal summer in the southern hemisphere. Other than that, one can check a myriad of cool things, and it is very fun to play around with the various graphs and results that can be achieved. For example, one can stop Earth’s rotation around itself, see what happens if it were tilted at 900, and compare the results for different latitudes.

One final note: I have not yet researched what the implications of having such a short precession period are. Axial precession is caused by torque, and makes Earth’s motion rather odd. The effects of such a motion on Earth’s inhabitants might be rather strong. It is also interesting to try and think what would cause such a quick movement, and if it can last for a large amount of time. Perhaps I will make another post of it one day.

The code, with default values simulating today’s Earth. For different values of parameters, just change the values given to the function “generateSolarFluxFunction”. All time units are in seconds, of course.

from math import pi, cos, sin, radians

def numericalIntegration(f, t0, t1, timestep):
    Numerically integrates the function f, returning the value
    of its definite integral from t0 to t1, using the trapezoidal rule.
    The parameter "timestep" determines the size of each trapezoid.
    @param f: A function of time.
    value1 = f(t0)
    currentTime = t0 + timestep
    total = 0
    while (currentTime <= t1):
        value2 = f(currentTime)

        total = total + timestep * (value1 + value2) * 0.5
        value1 = value2

        currentTime += timestep

    return total

def generateSolarFluxFunction(solarConstant = 1368.0,
        tiltDegrees = 23.3, latitudeDegrees = 0,
        dailyPeriod = 86400.0, yearlyPeriod = 365.25 * 86400.0,
        precessionPeriod = None):

    def solarFlux(t):
        latitude = radians(latitudeDegrees)
        tilt = radians(tiltDegrees)

        if dailyPeriod is not None:
            wd = 2.0 * pi / dailyPeriod
            wd = 0

        if yearlyPeriod is not None:
            wy = 2.0 * pi / yearlyPeriod
            wy = 0

        if precessionPeriod is not None:
            wp = 2.0 * pi / precessionPeriod
            wp = 0

        sinLat = sin(latitude)
        cosLat = cos(latitude)
        sinTilt = sin(tilt)
        cosTilt = cos(tilt)

        cosWd = cos(wd * t)
        sinWd = sin(wd * t)
        cosWp = cos(wp * t)
        sinWp = sin(wp * t)
        cosWy = cos(wy * t)
        sinWy = sin(wy * t)

        sunVector = [cosWy, sinWy, 0]

        patchVector = [
            (-cosLat * sinWd * cosTilt - sinLat * sinTilt) * cosWp
                - cosLat * cosWd * sinWp,
            cosLat * cosWd * cosWp +
                (-cosLat * sinWd * cosTilt - sinLat * sinTilt) * sinWp,
            sinLat * cosTilt - cosLat * sinWd * sinTilt,

        res = sum(sunVector[i] * patchVector[i] for i in xrange(3))
        if res < 0:
            return 0
        return solarConstant * res

    return solarFlux

def getDailyFluxResults():
    day = 86400.0
    time = 0

    # Change parameters here in order to get different effects.
    mySolarFlux = generateSolarFluxFunction()

    results = []
    while time < (366.0 * day):
        totalEnergyForThisDay = \
            numericalIntegration(mySolarFlux, time, time + day, 100.0)

        meanEnergyForThisDay = totalEnergyForThisDay / day        

        time += day

    return results


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s