Skip to content

2.2 Modeling the Seasonal Cycle

Learning Objectives

After completing this lecture, students should be able to:

  • Understand the concept of the seasonal cycle in surface air temperature
  • Interpret the parameters of a mathematical model representing the seasonal cycle
  • Estimate model parameters from observed temperature data
  • Run a Python script to simulate the seasonal temperature cycle
  • Analyze the simulation results and evaluate model performance using RMSE
  • Develop and evaluate a simple mathematical model of the observed seasonal temperature cycle.

Concept

As shown in the previous lecture, the observed surface air temperature in Utah exhibits a clear seasonal cycle. Because the seasonal cycle arises from the Earth's motion around the Sun, the temperature evolution over the year often resembles a sine or cosine curve.

In this section, we introduce a simple mathematical function to represent the seasonal cycle and compare it with the observed temperature data.

Using a mathematical representation allows us to describe the seasonal variability with a small number of parameters, which will be useful for later analyses in this course.

Model / Equation

The seasonal cycle of air temperature can be approximately represented by a sine or cosine function. If we denote time by \(t\), the annual temperature cycle \(T(t)\) can be written as

\[ T(t) = a_0 \sin\bigl(\omega (t - \phi)\bigr) + b_0 \]

or

\[ T(t) = a_0 \cos\bigl(\omega (t - \phi)\bigr) + b_0 \]

where

  • \(a_0\) : amplitude of the seasonal cycle
  • \(\omega\) : angular frequency
  • \(\phi\) : phase shift
  • \(b_0\) : mean temperature offset

The angular frequency determines the length of the cycle. For a seasonal cycle with a period of one year (365 days),

\[ \omega = \frac{2\pi}{365} \]

In this course, the time variable \(t\) is typically measured in days, so one full seasonal cycle corresponds to approximately 365 days.

Note that sine and cosine functions differ only by a phase shift of 90° (or \(\pi/2\) radians). Therefore, both equations represent the same seasonal cycle if the phase parameter \(\phi\) is adjusted accordingly.

Sample Code

Dataset (CSV files)

Surface Air Temperature (°F) at:

Original datasets can also be obtained from the Utah Climate Center:

https://climate.usu.edu/mchd/

Students may download additional weather station data from the Utah Climate Center and apply the same workflow to those datasets.

Python Script

We will use the template script developed in the previous section.

This script reads the weather station data and visualizes the temperature time series, which will serve as the observational reference when we compare it with the mathematical seasonal-cycle model.

The final code in this topic is here:

Exercise

In this exercise, we will construct a simple mathematical model of the seasonal cycle and gradually compare it with the observed temperature data.

Step 1. Define a Simple Time Coordinate

The dataset currently uses calendar date as the horizontal axis. For mathematical modeling, it is often more convenient to use a simple time index.

Create a new version of the script:

$ cp 2_1-read_temp_template.py 2_2-step1.py

Define a simple time coordinate in the Time setting section, just after the date definition:

t = xr.DataArray(range(len(fin)), dims=('time',), name='time', coords={'date': ('time', date)})

Here, t is the numerical time coordinate used for the mathematical model.
It takes the values:

  • 0, 1, 2, 3, ...

and represents time in days from the beginning of the dataset.

At the same time, we attach the original date information as an additional coordinate:

coords={'date': ('time', date)}

This is useful because:

  • time is convenient for mathematical expressions such as \(\cos(\omega t)\)
  • date preserves the original calendar information for interpretation and plotting

Then modify the temperature variable so that it uses the new time coordinate:

#temp = xr.DataArray(fin['temp'], dims=('date',), coords={'date': date})
temp = xr.DataArray(fin['temp'], dims=('time',), coords={'time': t, 'date': ('time', date)})

Here, temp is now defined on the time dimension rather than directly on the calendar-date dimension.
This allows us to:

  • use t for the mathematical model
  • keep date as a reference coordinate
  • maintain a connection between model time and calendar time

Update the plot so that the horizontal axis uses t instead of date:

#ax.plot(date, temp, 'o', markersize=2, color='C1')
ax.plot(t, temp, 'o', markersize=2, color='C1')

Run the script and confirm that the temperature time series is now plotted against a simple time index.

If you would like, you can also add grid lines to improve readability:

#ax.grid(True)
ax.grid(which='major', linewidth=1)
ax.grid(which='minor', linestyle='--', alpha=0.4)

This step is important because it separates:

  • calendar time (date) used to describe the observations
  • model time (t) used in the mathematical function

This distinction will be useful in the following steps when we construct the seasonal-cycle model.

Step 2. Add a Simple Sine Curve

Next, we introduce a simple mathematical function:

\[ y = a \sin(\omega (t - \phi)) + b \]

Start with simple parameters:

  • \(a = 1\)
  • \(\omega = 1\)
  • \(\phi = 0\)
  • \(b = 0\)

Create a new version of the script:

$ cp 2_2-step1.py 2_2-step2.py

Add the following import near the top of the script:

import numpy as np

Then add the following code before plotting:

w = 1.0
phi = 0.0

a = 1.0
b = 0.0

y = a * np.sin(w * (t - phi)) + b

Plot the curve together with the observations:

ax.plot(t, temp, 'o', markersize=2, color='C1', label='Observation')
ax.plot(t, y, '-', color='k', linewidth=2, label='Estimation')

This curve will not match the observations yet, but it demonstrates the behavior of the mathematical function.

Step 3. Construct a Seasonal Cycle Model

The seasonal cycle has a period of one year (365 days).
Therefore, the angular frequency of the model is

\[ \omega = \frac{2 \pi}{365} \]

Create a new version of the script:

$ cp 2_2-step2.py 2_2-step3.py

Then update the parameter in the script:

w = 2.0 * np.pi / 365.0

Next represent the seasonal cycle using a cosine function:

\[ y = a \cos(\omega (t - \phi)) + b \]

Modify the model equation:

y = a * np.cos(w * (t - phi)) + b

Now experiment with the parameters:

  • Amplitude (\(a\)) controls the seasonal temperature range
  • Offset (\(b\)) controls the mean temperature
  • Phase (\(\phi\)) controls the timing of the cycle (days)

Test different values of \(a\), \(b\), and \(\phi\), including negative and positive values, and observe how the curve changes. Run the script again and adjust these values until the curve roughly matches the observed seasonal cycle.

Step 4. Estimate the Parameters (Simple Approach)

There are several ways to estimate the model parameters. We will evaluate the performance of each estimate using the root-mean-square error (RMSE). A smaller RMSE indicates a better fit.

Create a new version of the script:

cp 2_2-step3.py 2_2-step4.py

The first estimate is a simple manual estimate. Surface air temperature typically reaches a minimum during mid-winter (around January 15) and a maximum during mid-summer (around July 15). Therefore, a reasonable first guess for the phase parameter is

phi = 15.0

Start by prescribing a and b as follows:

a = -20.0
b = 50.0

Run the script and compare the estimated seasonal cycle with the observations. Then adjust a, b, and phi to see how the model fit changes.

After defining the estimated curve y, include the RMSE calculation:

se = (temp - y) * (temp - y)
mse = se.mean(dim='time', skipna=True).item()
rmse = np.sqrt(mse)
print("RMSE =", rmse)

Or, equivalently, use the compact form:

rmse = np.sqrt(((temp - y)**2).mean(dim='time', skipna=True).item())
print("RMSE =", rmse)

The RMSE value can also be included in the figure title. Use standard Matplotlib commands for the title, axis labels, and legend:

title_name = "Seasonal Cycle Fit"

ax.set_xlabel('Time [days]')
ax.set_ylabel('Air Temperature [°F]')
ax.set_title(f'{title_name}\nRMSE = {rmse:.2f}°F')
ax.legend()

Alternatively, the RMSE can be shown as text inside the figure:

ax.set_xlabel('Time [days]')
ax.set_ylabel('Air Temperature [°F]')
ax.set_title(title_name)

ax.text(
    0.98, 0.95,
    f'RMSE = {rmse:.2f}°F',
    transform=ax.transAxes,
    ha='right',
    va='top'
)

ax.legend()

Run the code and check the RMSE output.

Another simple method to estimate the parameters a and b is to use the mid-summer and mid-winter temperatures. Surface air temperature typically reaches a maximum during mid-summer (around July 15) and a minimum during mid-winter (around January 15). Replace the parameter estimation with the following:

t_mod = t % 365

temp_summer = temp.where((t_mod >= 190) & (t_mod <= 200), drop=True)
temp_winter = temp.where((t_mod >= 10) & (t_mod <= 20), drop=True)

T_summer = temp_summer.mean(dim='time', skipna=True).item()
T_winter = temp_winter.mean(dim='time', skipna=True).item()

print("Mid-summer temperature =", T_summer)
print("Mid-winter temperature =", T_winter)

b = 0.5 * (T_summer + T_winter)
a = -0.5 * (T_summer - T_winter)

Run the code again and compare the RMSE between the two estimation methods.

Step 5. Estimate the Parameters (Least Squares Method)

One of the most common methods for estimating model parameters is the least-squares method. This method assumes that the observed values can be approximated by the model values. In the present example, we assume that

\[ y_{\mathrm{obs}} \approx a \cos(\omega (t - \phi)) + b \]

To simplify the notation, define

\[ f(t) = \cos(\omega (t - \phi)) \]

Then the model can be written as

\[ y_{\mathrm{obs}} \approx a f(t) + b \]

If we move \(b\) to the left-hand side and multiply both sides by \(f(t)\), we obtain

\[ \bigl(y_{\mathrm{obs}} - b\bigr) f(t) = a\, f(t)\, f(t) \]

Taking the average of both sides gives

\[ a = \frac{\left< \bigl(y_{\mathrm{obs}} - b\bigr) f(t) \right>} {\left< f(t)\, f(t) \right>} \]

where \(\langle \cdot \rangle\) denotes the average over time.
The parameter \(b\) can be estimated as the mean of \(y_{\mathrm{obs}}\).

Create a new version of the script:

$ cp 2_2-step4.py 2_2-step5.py

The Python code corresponding to this formula is:

temp_mean = temp.mean(dim='time', skipna=True).item()
b = temp_mean
print("Mean temperature, b =", b)

f = np.cos(w * (t - phi))
a = ((temp - b) * f).mean(dim='time', skipna=True).item() / (f * f).mean(dim='time', skipna=True).item()

print("a =", a)

After estimating \(a\) and \(b\), run the code again and compare the RMSE with the previous estimation methods.

Discussion

After completing the steps above, consider the following questions:

  • Which parameters (\(a\), \(\omega\), \(\phi\), \(b\)) control the amplitude, frequency, phase, and mean state of the seasonal cycle?
  • How does the phase parameter \(\phi\) shift the seasonal cycle forward or backward in time?
  • How well can a simple sine or cosine function represent the observed temperature variability?
  • What are the limitations of this simple model in describing real atmospheric temperature variations?

This model provides a first-order mathematical representation of the seasonal cycle, forming the foundation for more advanced fitting methods such as the least-squares approach.

Key Points

  • The seasonal cycle of temperature can be approximated using a sine or cosine function.
  • The parameters have clear physical interpretations:
  • \(a\): amplitude (seasonal range)
  • \(\omega\): frequency (annual cycle)
  • \(\phi\): phase (timing of peaks and troughs)
  • \(b\): mean temperature (background state)
  • Using a simple mathematical model allows us to quantify and interpret seasonal variability.
  • Model performance can be evaluated using RMSE, providing an objective measure of fit.

Homework

Complete the following tasks:

  1. Apply the seasonal cycle model to the Venice temperature dataset (Venice_temp.csv).
  2. Estimate the parameters (\(a\), \(\phi\), \(b\)) using any method (e.g., manual tuning, mean value, seasonal extremes, or least-squares method).
  3. Generate a figure that includes: - Observed temperature
    - Estimated seasonal cycle
    - RMSE displayed on the figure
  4. Briefly describe: - How you determined the parameters
    - How well the model fits the observations based on RMSE and visual comparison