Skip to content

2.3 Soil Temperature and External Forcing

Learning Objectives

After completing this lecture, students should be able to:

  • Understand the decomposition of observed temperature into climatology, external forcing, and internal variability.
  • Interpret the physical meaning of the seasonal cycle, linear trend, and residual variability.
  • Estimate seasonal and linear-trend parameters using least-squares methods and polynomial fitting.
  • Run Python scripts to process observational data and simulate soil-temperature variability.
  • Analyze soil-temperature observations by removing the seasonal cycle and trend and examining internal variability across different time scales.

Concept

Previously, we decomposed observed temperature into a climatological mean and an anomaly:

\[ T = \overline{T} + T' \]

where:

  • \(T\): observed temperature
  • \(\overline{T}\): climatological temperature
  • \(T'\): temperature anomaly

Temperature anomalies can be further separated into variability associated with external forcing and internal variability:

\[ T' = T_{\mathrm{ext}} + T_{\mathrm{int}} \]

Therefore,

\[ \boxed{ T = T_{\mathrm{clm}} + T_{\mathrm{ext}} + T_{\mathrm{int}} } \]

where:

  • \(T_{\mathrm{clm}}\): climatological seasonal cycle
  • \(T_{\mathrm{ext}}\): externally forced change, approximated here by a long-term trend
  • \(T_{\mathrm{int}}\): internal variability remaining after removing the climatology and trend

This decomposition provides a simple framework for separating different sources of observed temperature variability.

Model / Equation

We approximate the climatological seasonal cycle using the cosine function introduced previously:

\[ T_{\mathrm{clm}}(t) = a \cos[\omega(t-\phi)]+b \]

We then represent the externally forced component using a linear trend:

\[ T_{\mathrm{ext}}(t)=ct \]

Combining these components gives:

\[ T(t) = a \cos[\omega(t-\phi)]+b +ct +T_{\mathrm{int}}(t) \]

or

\[ \boxed{ T(t)=T_{\mathrm{clm}}(t)+T_{\mathrm{ext}}(t)+T_{\mathrm{int}}(t) } \]

The internal variability can therefore be estimated as the residual:

\[ \boxed{ T_{\mathrm{int}}(t) = T(t)-T_{\mathrm{clm}}(t)-T_{\mathrm{ext}}(t) } \]

where:

  • \(a\): amplitude of the seasonal cycle
  • \(\omega\): angular frequency of the annual cycle
  • \(\phi\): phase, representing the timing of the seasonal cycle
  • \(b\): mean background temperature
  • \(c\): linear trend coefficient
  • \(t\): time
  • \(T_{\mathrm{int}}(t)\): residual/internal variability

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 formulate the mathmatical function, simulating seasonal cycle and linear trend of soil temperature.

Step 1. Read an Advanced CSV File

Download the station observation dataset DugwayUT.csv, which contains multiple variables, including air temperature and soil temperature.

Create a new version of the script:

$ cp 2_2-mathfunc_template.py 2_3-soil_step1.py

Modify the CSV file path and title:

#file_csv = 'greenfarm_temp.csv'
#title_name = 'Greenville Farm Weather Station'
#file_csv = 'Venice_temp.csv'
#title_name = 'Venice Weather Station'
file_csv = 'DugwayUT.csv'
title_name = 'Dugway, UT'

Run the code. You will likely encounter an error when reading the CSV file. This is expected because the file structure is more complex than the previous datasets.

Inspect the CSV File:

To understand the file structure, inspect the CSV file directly from the command line. A CSV file is a text-based file, and you can view its contents using the less command. The -N option displays line numbers:

$ less -N DugwayUT.csv

Useful commands within less:

  • Enter → move forward by one line
  • Space → move forward by one page
  • b → move backward by one page
  • ↑ / ↓ → move one line up or down
  • q → quit

By inspecting the file, you will notice that:

  • The file contains metadata in the first several lines
  • The actual data table begins around line 72
  • The dataset includes many columns

Update the CSV Reading Method:

To correctly read the file, modify the read_csv function to skip the metadata and use the appropriate header line:

fin = pd.read_csv(
    file_csv,
    na_values=['M', 'T'],
    encoding='ISO-8859-1',
    delimiter=',',
    header=71,
    #names=['date', 'temp']
)

Note that the line numbers displayed by the less command start from 1, whereas Python (and pandas) uses zero-based indexing.

Therefore, if the header appears on line 72 in less, the corresponding index in pandas is:

\[ 72 - 1 = 71 \]

This is why we set:

header=71

Adjust Column Names:

Run the code again. You may encounter another error indicating that the column names are not recognized. This happens because the variable names differ from the simpler datasets used previously.

Update the variable definitions using the correct column names:

fin['date'] = pd.to_datetime(fin['Date'])
fin['temp'] = pd.to_numeric(fin['Air Temperature Maximum (degF)'])

Final Check:

Run the code again. If everything is set correctly, the script should execute without errors and produce the temperature figure.

This step demonstrates how to handle real-world observational datasets, which often include:

  • metadata headers
  • multiple variables
  • non-standard formatting

Understanding how to inspect and preprocess such data is an essential skill in climate data analysis.

Step 2. Analyze Soil Temperature

Next, examine the soil temperature and characterize its temporal evolution.

Create a new version of the script:

$ cp 2_3-soil_step1.py 2_3-soil_step2.py

Then, read the soil temperature at a depth of 4 inches and update the y-axis label:

fin['temp'] = pd.to_numeric(
    fin['Soil Temperature Observed -4in (degF) Mean of Hourly Values']
)

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

Run the code and examine the figure. Compare the temporal evolution of soil temperature with that of air temperature.
Are there differences in the timing and magnitude of maximum and minimum temperatures?

Remove the seasonal cycle using a running mean:

Because the seasonal cycle dominates the signal, it can obscure lower-frequency variability.
One simple way to reduce the seasonal cycle is to apply a 365-day running mean.

Add the following calculation before plotting, and include it in the figure:

temp_rm = temp.rolling(time=365, center=True).mean()

ax.plot(t, temp_rm, '-', linewidth=2, color='C0', label='Running mean')
#ax.plot(t, y, '-', color='k', linewidth=2, label='Estimation')

This smoothed time series suppresses the seasonal cycle and helps reveal longer-term variations or trends.

Remove the seasonal cycle using anomalies:

In climate science, variability is often analyzed using anomalies, defined as deviations from a climatological mean.

Since we previously estimated the seasonal cycle using a cosine function, we can define temperature anomalies by removing this estimated climatology:

#temp_rm = temp.rolling(time=365, center=True).mean()
temp_anm = temp - y
temp_anm_rm = temp_anm.rolling(time=365, center=True).mean()

#ax.plot(t, temp, 'o', markersize=2, color='C1', label='Observation')
#ax.plot(t, temp_rm, '-', linewidth=2, color='C0', label='Running mean')
#ax.plot(t, y, '-', color='k', linewidth=2, label='Estimation')
ax.plot(t, temp_anm, 'o', markersize=2, color='C1', label='Anomaly')
ax.plot(t, temp_anm_rm, '-', linewidth=2, color='C0', label='Running mean')

Run the code and examine how the temperature variability appears after removing the seasonal cycle.

Finally, experiment with different running-mean windows (e.g., time=30, 90, 120, 180) and discuss how the apparent variability changes with the smoothing scale.

Step 3. Estimate the Linear Trend

In this step, we estimate the linear trend in soil temperature.

We represent the trend component as

\[ T_{\mathrm{trend}}(t) = c_0 t + d_0 \]

However, in the current model we have already estimated the mean offset term \(b\). Therefore, we first focus on estimating only the slope term and write the model as

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

where:

  • \(a\) : amplitude of the seasonal cycle
  • \(\omega\) : angular frequency
  • \(\phi\) : phase
  • \(b\) : mean offset
  • \(c\) : linear trend coefficient

Create a new version of the script:

$ cp 2_3-soil_step2.py 2_3-soil_step3.py

Estimate the trend visually first:

Run the copied script with a 365-day running mean:

temp_anm_rm = temp_anm.rolling(time=365, center=True).mean()

From the graph, estimate the slope of the blue line. For example, you may observe that the temperature anomalies increase by about 8°F over 5 years.

Before using this estimate, first modify the code to include the linear trend term. Set c = 0.0, then update the model equation. We also need to restore the plot for the observation and estimation:

c = 0.0
print("c =", c)

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

#ax.plot(date, temp, 'o', markersize=2, color='C1')
ax.plot(t, temp, 'o', markersize=2, color='C1', label='Observation')
#ax.plot(t, temp_anm, 'o', markersize=2, color='C1', label='Observation')
#ax.plot(t, temp_rm, '-', linewidth=2, color='C0', label='Running mean')
#ax.plot(t, temp_anm_rm, '-', linewidth=2, color='C0', label='Running mean')
ax.plot(t, y, '-', color='k', linewidth=2, label='Estimation')

Run the code and take note of the RMSE value.

Next, set c using your visual estimate. You also need to adjust the temperature residual used to estimate the seasonal amplitude by subtracting both the mean and the linear trend:

c = 8.0 / (5.0 * 365.0)

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

Run the code again and compare the RMSE value. You may find that the RMSE becomes larger, indicating that the visual estimate of the trend is not optimal.

Estimate the slope using least squares:

We can apply the least-squares method to estimate the trend coefficient.

Previously, we used least squares to estimate the seasonal amplitude. Here we use the same idea, but replace the function with \(f_1(t) = t\) to estimate the slope \(c\).

Add the following lines:

f1 = t
c = ((temp - b) * f1).mean(dim='time', skipna=True).item() / (f1 * f1).mean(dim='time', skipna=True).item()
print("c =", c)

Run the code and assess the RMSE value again.

We still need a consistent estimate of the intercept. Introduce the mean of time and temperature, and revise the estimation:

t_mean = t.mean().item()
temp_mean = temp.mean(dim='time', skipna=True).item()

f1 = t
c = ((temp - temp_mean) * f1).mean(dim='time', skipna=True).item() / (f1 * f1).mean(dim='time', skipna=True).item()
print("c =", c)

b = temp_mean - c * t_mean
print("Intersection, b =", b)

Run the code and assess the RMSE value again.

Estimate the slope using a linear fitting function (polynomial fit):

Another common and convenient approach is to use a built-in linear fitting method. In Python, this can be done using numpy.polyfit, which directly fits a first-order polynomial:

\[ T(t) = c t + d \]

Add the following lines:

#f1 = t
#c = ((temp - temp_mean) * f1).mean(dim='time', skipna=True).item() / (f1 * f1).mean(dim='time', skipna=True).item()
#print("c =", c)
#b = temp_mean - c * t_mean
#print("Intersection, b =", b)

coef = np.polyfit(t.values, temp.values, 1)
c = coef[0]
b = coef[1]

print("c =", c)
print("Intersection, b =", b)

Run the code and compare the result with the least-squares estimate.

  • Are the slope values similar?
  • Does the fitted line better represent the long-term behavior?
  • How does the RMSE compare?

Difference between least-squares method and polynomial fit:

Although both approaches rely on the same least-squares principle, their implementations differ slightly.

  • The manual least-squares method explicitly applies the mathematical formula using averages (covariance and variance). This helps clarify how the slope and intercept are derived from the data.

  • The np.polyfit function performs the same least-squares optimization internally, but in a more automated and general way. It estimates both slope and intercept simultaneously and can be extended to higher-order polynomial fits.

When the least-squares method is implemented consistently (including proper estimation of both slope and intercept), both approaches should produce nearly identical results.

In practice: - Use the manual method to understand the mathematics
- Use np.polyfit for efficiency and flexibility in real applications

Step 4. Estimate the Internal Variability

Typically, daily temperature contains a seasonal cycle and a long-term linear trend. These two components are common sources of external forcing in the climate system. Once we estimate and remove these components, we can better assess the internal variability.

Create a new version of the script:

$ cp 2_3-soil_step3.py 2_3-soil_step4.py

Compare how the anomalies look when we remove:

  • only the seasonal cycle
  • both the seasonal cycle and the linear trend

Here, y1 includes only the seasonal cycle, whereas y2 includes both the seasonal cycle and the linear trend. The anomalies are defined as deviations from these estimates.

We do not need RMSE in this step, so comment out those lines and remove the right title.

y1 = a * np.cos(w * (t - phi)) + b
y2 = a * np.cos(w * (t - phi)) + b + c * t
y_trend = c * t

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

temp_anm1 = temp - y1
temp_anm2 = temp - y2

ax.plot(t, temp_anm1, 'o', markersize=2, color='C1', label='No seasonal cycle')
ax.plot(t, y_trend, '-', color='k', linewidth=2, label='Linear trend')
ax.plot(t, temp_anm2, 'o', markersize=2, color='C0', label='No seasonal cycle & detrend')

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

Run the code and discuss the outcome.

Also try both the least squared method and polynominal fit approach, then observe the outcomes.

Visualize internal variability:

By applying a running mean, internal variability becomes more visible. Apply a running mean to the anomalies and show the result:

temp_anm1_rm = temp_anm1.rolling(time=365, center=True).mean()
temp_anm2_rm = temp_anm2.rolling(time=365, center=True).mean()

#ax.plot(t, y_trend, '-', color='k', linewidth=2, label='Linear trend')
ax.plot(t, temp_anm1_rm, '-', linewidth=2, color='C1', label='Running mean 1')
ax.plot(t, temp_anm2_rm, '-', linewidth=2, color='C0', label='Running mean 2')

This figure may still look busy. To improve the visualization, use lighter colors for the points and darker colors for the lines:

ax.plot(t, temp_anm1, 'o', markersize=2, color='pink', label='No seasonal cycle')
ax.plot(t, temp_anm2, 'o', markersize=2, color='lightblue', label='No seasonal cycle & detrend')
#ax.plot(t, y_trend, '-', color='k', linewidth=2, label='Linear trend')
ax.plot(t, temp_anm1_rm, '-', linewidth=2, color='red', label='Running mean 1')
ax.plot(t, temp_anm2_rm, '-', linewidth=2, color='darkblue', label='Running mean 2')

Finally, apply several running-mean windows, such as:

  • submonthly (time=15)
  • monthly (time=31)
  • seasonal (time=61)
  • interannual (time=91)

and discuss how the apparent variability changes with the smoothing scale.

Key Points

  • Observed temperature can be decomposed into three main components:
  • Seasonal cycle (periodic variation)
  • Long-term trend (linear change over time)

  • Internal variability (residual fluctuations)

  • A simple mathematical model combining these components is:

$$ T(t) = a \cos(\omega (t - \phi)) + b + c t $$

  • The parameters have clear physical interpretations:
  • \(a\) : amplitude of the seasonal cycle
  • \(\omega\) : frequency (annual cycle)
  • \(\phi\) : phase (timing of peaks and troughs)
  • \(b\) : mean background temperature

  • \(c\) : long-term trend

  • Removing the seasonal cycle and trend allows us to isolate internal variability, which is essential for understanding climate fluctuations.

  • Running mean smoothing helps reveal low-frequency signals, but the choice of window length affects the interpretation.

Homework

Complete the following tasks:

  1. Analyze soil temperature at different depths

Use the Dugway dataset to extract soil temperature at:

  • 8 inches
  • 20 inches
  • 40 inches
  1. Compare key components

For each depth, analyze and compare:

  • Seasonal cycle (amplitude and phase)
  • Linear trend (magnitude and sign)
  • Internal variability (after removing seasonal cycle and trend)
  1. Interpret the results

Discuss how soil depth influences:

  • the strength of the seasonal cycle
  • the timing (phase delay)
  • the magnitude of variability
  1. Method comparison

Estimate the linear trend using:

  • the least-squares method (as in the exercise)
  • a direct linear fitting method (e.g., numpy.polyfit)

Compare the results and discuss: - differences in estimated slope
- differences in RMSE
- which method provides a better fit and why

  1. Submit
  • At least one figure for each depth showing:
    • observation
    • fitted seasonal cycle
    • linear trend
    • anomalies (optional but encouraged)
  • A short written discussion summarizing your findings