Skip to content

7.2 External Forcing I: Seasonal Cycle

Learning Objectives

After completing this lecture, students should be able to:

  • distinguish external forcing from internal variability,
  • represent the seasonal cycle using a periodic forcing,
  • understand how periodic forcing modifies the climate state,
  • perform climate model experiments by modifying only the external forcing, and
  • interpret the combined effects of internal variability and seasonal forcing.

Concept

In the previous section, the Lorenz climate model simulated internal variability without any external forcing. Real climate systems, however, are continually influenced by changes in the Earth's energy balance. One of the most important periodic forcings is the seasonal cycle of net surface radiation, which varies as the Earth revolves around the Sun.

To represent this process, we introduce an externally prescribed forcing into the climate model. Unlike the internal dynamics, the external forcing depends only on time and is independent of the model state.

Although the Lorenz model is highly idealized, this simple forcing provides a useful framework for understanding how periodic external forcing interacts with internally generated climate variability.

Model / Equation

As in the previous section, the total climate state is decomposed into internal and externally forced components,

\[ X_{\rm total} = X_{\rm int} + X_{\rm ext}. \]

The governing equation becomes

\[ \frac{dX}{dt} = F_{\rm int}(X) + F_{\rm ext}(t), \]

where

  • \(F_{\rm int}\) represents the nonlinear Lorenz dynamics,
  • \(F_{\rm ext}\) represents an externally prescribed forcing.

To represent the seasonal cycle, the external forcing is approximated by a cosine function,

\[ F_x(t) = A_x \cos \left( \frac{2\pi t}{P} \right), \]

where

  • \(A_x\) is the forcing amplitude,
  • \(P\) is the forcing period.

Similarly,

\[ F_y(t) = A_y \cos \left( \frac{2\pi t}{P} \right), \]
\[ F_z(t) = A_z \cos \left( \frac{2\pi t}{P} \right). \]

The periodic forcing is implemented in

def external_forcing_tendency(TIME):

For example,

PERIOD = 12.0
OMEGA = 2.0 * np.pi / PERIOD

Fx = AMP_X * np.cos(OMEGA * TIME)
Fy = AMP_Y * np.cos(OMEGA * TIME)
Fz = AMP_Z * np.cos(OMEGA * TIME)

return Fx, Fy, Fz

The forcing amplitudes are prescribed in run_model_climate-v3.sh, for example,

AMP_X="-20.0"
AMP_Y="0.0"
AMP_Z="0.0"

Different climate experiments are performed simply by modifying the forcing amplitudes while leaving the model equations unchanged.

Sample Code

The source code for this lecture is provided as the compressed archive

Download the archive and extract it using

tar -xzvf Model_src-v3.0.tar.gz

This creates the directory

Model_src-v3/
├── figures/
│   └── fig-timeseries_climate.py
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

The individual source files are also available below.

Compared with the previous lecture, only the external forcing has been modified. The model structure, numerical integration, and visualization scripts remain unchanged.

Exercise

Step 1. Environment Setting

Create a new experiment directory for this section,

mkdir Model_EXP7_2

Move into the experiment directory,

cd Model_EXP7_2

Copy the model source files,

cp ../Model_src-v3/*.* .

Verify that the directory contains

Model_EXP7_2/
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

Step 2. Model Spin-Up

Open run_model_climate-v3.sh and configure a control simulation without external forcing,

DEBUG_OUTPUT="False"

REST="False"

TSTR="0"
TEND="60"

AMP_X="0.0"
AMP_Y="0.0"
AMP_Z="0.0"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

After the simulation finishes, verify that the output and restart files have been created,

ls out

ls rest

Inspect one of the model output files,

ncdump out/out.0000.nc

Questions

  • What variables are stored in the output file?
  • Why is a model spin-up necessary?
  • Why are restart files created during the simulation?

Step 3. Seasonal Forcing Run

Next, continue the simulation from the restart file and apply a seasonal forcing.

Open run_model_climate-v3.sh and modify

DEBUG_OUTPUT="False"

REST="True"

TSTR="60"
TEND="180"

AMP_X="-20.0"
AMP_Y="0.0"
AMP_Z="0.0"

Then modify the external forcing function in Lmodel_climate-v3.py:

def external_forcing_tendency(TIME):

    PERIOD = 12.0
    OMEGA = 2.0 * np.pi / PERIOD

    Fx = AMP_X * np.cos(OMEGA * TIME)
    Fy = AMP_Y * np.cos(OMEGA * TIME)
    Fz = AMP_Z * np.cos(OMEGA * TIME)

    return Fx, Fy, Fz

Run the model again,

./run_model_climate-v3.sh

Inspect the output files,

ls out

Questions

  • Which Lorenz variable is directly forced?
  • Why is the restart file used instead of restarting from the initial condition?
  • What is the period of the prescribed forcing?

Step 4. Visualization

Create a directory for visualization scripts,

mkdir figures

Copy the time-series plotting script,

cp ../Model_src-v3/figures/fig-timeseries_climate.py figures/

Move to the figures directory,

cd figures

Generate the time-series figure,

python fig-timeseries_climate.py

If necessary, modify the plotting period to display the complete simulation,

TSTR = 0
TEND = 180

Run the plotting script again.

Compare the time series before and after the seasonal forcing is applied at t = 60.

Pay particular attention to the differences among the three Lorenz variables.

Questions

  • How does the behavior of X change after the seasonal forcing is introduced?
  • Are similar periodic variations immediately visible in Y and Z?
  • Why is the seasonal signal more apparent in X than in the other variables?
  • Does the seasonal forcing eliminate the chaotic fluctuations?
  • From the total climate state alone, is it possible to distinguish the externally forced response from the internally generated variability?

Step 5. Long-Term Control Simulation

To better identify the seasonal cycle, perform a longer control simulation using the restart file generated during the spin-up experiment.

Create a new experiment directory,

mkdir CTL

Copy the model source files,

cp *.py *.sh CTL/

Move into the control experiment directory,

cd CTL

Create a directory for the restart files,

mkdir rest

Copy the restart file from the previous experiment,

cp ../rest/restart.0060.nc rest/

Open run_model_climate-v3.sh and configure the experiment as

DEBUG_OUTPUT="True"

REST="True"

TSTR="60"
TEND="240"

AMP_X="-20.0"
AMP_Y="-15.0"
AMP_Z="-25.0"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

Create a directory for the visualization scripts,

mkdir figures

Copy the time-series plotting script,

cp ../../Model_src-v3/figures/fig-timeseries_climate.py figures/

Move to the figures directory,

cd figures

Modify the plotting period,

TSTR = 60
TEND = 240

Generate the time-series figure,

python fig-timeseries_climate.py

Step 6. Estimating the Seasonal Cycle

The total climate state consists of both the externally forced seasonal cycle and internally generated variability. To estimate the externally forced response, we first rearrange the model output into a Year × Month format.

Each row represents one model year, and each column represents the same month within the seasonal cycle. The climatological seasonal cycle is obtained by averaging across the year dimension.

Step 6.1 Plot Individual Seasonal Cycles

Copy the seasonal-cycle plotting script,

cp ../../Model_src-v3/figures/fig-seasonal_cycle-template.py ./fig-seasonal_cycle-raw.py

Open the script and confirm:

TSTR = 60
TEND = 240

PERIOD = 12.0
N_YEAR = 3

Run the script,

python fig-seasonal_cycle-raw.py

Questions:

  • Are the three seasonal cycles identical?
  • Why do individual years differ even though the same seasonal forcing is applied every year?
  • Which variable shows the clearest seasonal cycle?

Step 6.2 Estimate the Climatological Seasonal Cycle

Create a new plotting script,

cp fig-seasonal_cycle-template.py fig-seasonal_cycle-clm.py

Open fig-seasonal_cycle_clm.py and compute the climatological seasonal cycle by averaging over the year dimension. Put the following code just before the figure plot,

Xclim = np.mean(Xyear, axis=0)
Yclim = np.mean(Yyear, axis=0)
Zclim = np.mean(Zyear, axis=0)

Modify the figure so that individual years are plotted with thin, light-gray lines,

for iyear in range(N_YEAR):

    ax[0].plot(month, Xyear[iyear, :], color="lightgray", linewidth=0.8)
    ax[1].plot(month, Yyear[iyear, :], color="lightgray", linewidth=0.8)
    ax[2].plot(month, Zyear[iyear, :], color="lightgray", linewidth=0.8)

Then overlay the climatological seasonal cycle using a thick line,

ax[0].plot(month, Xclim, color="C0", linewidth=3.0, label="Climatology")
ax[1].plot(month, Yclim, color="C0", linewidth=3.0, label="Climatology")
ax[2].plot(month, Zclim, color="C0", linewidth=3.0, label="Climatology")

Run the script,

python fig-seasonal_cycle-clm.py

After confirming that the script works, increase the number of years,

N_YEAR = 15

and run the script again.

Questions:

  • How does the climatology change when N_YEAR increases from 3 to 15?
  • Why does averaging more years produce a smoother seasonal cycle?
  • Which variable exhibits the strongest climatological seasonal response?

Step 6.3 Evaluate the Estimated Seasonal Cycle

In this experiment, the model was run with

DEBUG_OUTPUT="True"

Therefore, the model output includes the true externally forced components:

  • Xext
  • Yext
  • Zext

These variables are generally unavailable in observations, but they provide an ideal reference for evaluating the climatological seasonal cycle estimated from the total climate state.

Create a new plotting script,

cp fig-seasonal_cycle-clm.py fig-seasonal_cycle-rmse.py

Read the external components,

Xext = fin.Xext.sel(time=slice(TSTR, TEND))
Yext = fin.Yext.sel(time=slice(TSTR, TEND))
Zext = fin.Zext.sel(time=slice(TSTR, TEND))

The externally forced components contain only the seasonal perturbation, whereas the total climate state also includes the mean background state. Before comparing the two, adjust the mean of each externally forced component so that it matches the mean of the corresponding total variable:

Xext = Xext + (X.mean() - Xext.mean())
Yext = Yext + (Y.mean() - Yext.mean())
Zext = Zext + (Z.mean() - Zext.mean())

This bias adjustment allows the comparison to focus on the seasonal cycle itself rather than differences in the mean state.

Next, rearrange Xext, Yext, and Zext into Year × Month arrays using the same procedure as for Xyear, Yyear, and Zyear.

Compute the climatological seasonal cycles,

Xext_clm = np.mean(Xext_year, axis=0)
Yext_clm = np.mean(Yext_year, axis=0)
Zext_clm = np.mean(Zext_year, axis=0)

Calculate the RMSE between the estimated climatology and the true externally forced climatology,

rmseX = np.sqrt(np.mean((Xclm - Xext_clm)**2))
rmseY = np.sqrt(np.mean((Yclm - Yext_clm)**2))
rmseZ = np.sqrt(np.mean((Zclm - Zext_clm)**2))

Overlay the true externally forced climatology using a dashed line,

ax[0].plot(month, Xext_clm, color="k", linestyle="--", linewidth=2.0,
           label="True external")

ax[1].plot(month, Yext_clm, color="k", linestyle="--", linewidth=2.0,
           label="True external")

ax[2].plot(month, Zext_clm, color="k", linestyle="--", linewidth=2.0,
           label="True external")

Display the number of years and the RMSE in each panel,

ax[0].text(
    0.02, 0.88,
    f"N_YEAR = {N_YEAR}\nRMSE = {rmseX:.2f}",
    transform=ax[0].transAxes,
    bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)

Repeat the analysis with different values of

N_YEAR

for example,

N_YEAR = 3

and

N_YEAR = 15

Run the script,

python fig-seasonal_cycle-rmse.py

Questions

  • How does the RMSE change as N_YEAR increases?
  • Why does averaging more years improve the estimate of the externally forced seasonal cycle?
  • Which variable has the largest RMSE?
  • Why is it important to adjust the mean before comparing the estimated and true seasonal cycles?

Step 7. Estimating Internal Variability

After estimating the climatological seasonal cycle, the internally generated variability can be estimated by removing the climatology from the total climate state:

\[ X' = X - X_{\rm clm}. \]

This is the same approach commonly used in climate-data analysis, where anomalies are computed by subtracting the climatological seasonal cycle.

Create a new plotting script from the time-series script,

cp fig-timeseries_climate.py fig-timeseries_internal_variability.py

In fig-timeseries_internal_variability.py, prescribe the analysis period and number of years:

import numpy as np

TSTR = 60
TEND = 240

PERIOD = 12.0
N_YEAR = 3

Add a check to ensure that the requested number of years is available:

n_available = int((TEND - TSTR) / PERIOD)

if N_YEAR > n_available:
    raise ValueError(
        f"N_YEAR ({N_YEAR}) exceeds the available number of years ({n_available}). "
        "Please reduce N_YEAR or increase TEND."
    )

Read both the total climate state and the true internal component from the debug output:

X = fin.X.sel(time=slice(TSTR, TEND))
Y = fin.Y.sel(time=slice(TSTR, TEND))
Z = fin.Z.sel(time=slice(TSTR, TEND))

Xint = fin.Xint.sel(time=slice(TSTR, TEND))
Yint = fin.Yint.sel(time=slice(TSTR, TEND))
Zint = fin.Zint.sel(time=slice(TSTR, TEND))

Then rearrange the total climate state into Year × Month arrays as in Step 6,

DT = float(time[1] - time[0])
N_MONTH = int(PERIOD / DT)
month = np.arange(N_MONTH) * DT

Xyear = np.zeros((N_YEAR, N_MONTH))
Yyear = np.zeros((N_YEAR, N_MONTH))
Zyear = np.zeros((N_YEAR, N_MONTH))

for iyear in range(N_YEAR):

    t0 = TSTR + iyear * PERIOD
    t1 = t0 + PERIOD

    Xplot = X.sel(time=slice(t0, t1 - 0.5 * DT))
    Yplot = Y.sel(time=slice(t0, t1 - 0.5 * DT))
    Zplot = Z.sel(time=slice(t0, t1 - 0.5 * DT))

    Xyear[iyear, :] = Xplot.values
    Yyear[iyear, :] = Yplot.values
    Zyear[iyear, :] = Zplot.values

Rearrange the total climate state into Year × Month arrays, estimate the climatological seasonal cycle, and remove it:

Xclm = np.mean(Xyear, axis=0)
Yclm = np.mean(Yyear, axis=0)
Zclm = np.mean(Zyear, axis=0)

Xanom = Xyear - Xclm[None, :]
Yanom = Yyear - Yclm[None, :]
Zanom = Zyear - Zclm[None, :]

Convert the anomalies back to one-dimensional time series:

Xanom_ts = Xanom.reshape(N_YEAR * N_MONTH)
Yanom_ts = Yanom.reshape(N_YEAR * N_MONTH)
Zanom_ts = Zanom.reshape(N_YEAR * N_MONTH)

Create the corresponding time coordinate,

time_anom = time.values[:N_YEAR * N_MONTH]

Prepare the true internal variability from the debug output using the same length:

Xint_ts = Xint.values[:N_YEAR * N_MONTH]
Yint_ts = Yint.values[:N_YEAR * N_MONTH]
Zint_ts = Zint.values[:N_YEAR * N_MONTH]

Xint_ts = Xint_ts - np.mean(Xint_ts)
Yint_ts = Yint_ts - np.mean(Yint_ts)
Zint_ts = Zint_ts - np.mean(Zint_ts)

Compute RMSE between the estimated and true internal variability:

rmseX = np.sqrt(np.mean((Xanom_ts - Xint_ts)**2))
rmseY = np.sqrt(np.mean((Yanom_ts - Yint_ts)**2))
rmseZ = np.sqrt(np.mean((Zanom_ts - Zint_ts)**2))

Open fig-timeseries_internal_variability.py.

Replace the original plotting commands,

ax[0].plot(time, X, color="C0", linewidth=1.0)
ax[1].plot(time, Y, color="C0", linewidth=1.0)
ax[2].plot(time, Z, color="C0", linewidth=1.0)

with

ax[0].plot( time_anom, Xanom_ts, color="C0", linewidth=1.0, label="Estimated")
ax[0].plot( time_anom, Xint_ts, color="k", linestyle="--", linewidth=1.0, label="True State")

ax[1].plot( time_anom, Yanom_ts, color="C0", linewidth=1.0, label="Estimated")
ax[1].plot( time_anom, Yint_ts, color="k", linestyle="--", linewidth=1.0, label="True State")

ax[2].plot( time_anom, Zanom_ts, color="C0", linewidth=1.0, label="Estimated")
ax[2].plot( time_anom, Zint_ts, color="k", linestyle="--", linewidth=1.0, label="True State")

Then add the RMSE and the number of years to each panel:

ax[0].text(
    0.02, 0.88,
    f"N_YEAR = {N_YEAR}\nRMSE = {rmseX:.2f}",
    transform=ax[0].transAxes,
    bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)

ax[1].text(
    0.02, 0.88,
    f"N_YEAR = {N_YEAR}\nRMSE = {rmseY:.2f}",
    transform=ax[1].transAxes,
    bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)

ax[2].text(
    0.02, 0.88,
    f"N_YEAR = {N_YEAR}\nRMSE = {rmseZ:.2f}",
    transform=ax[2].transAxes,
    bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)

The resulting figure compares the estimated internal variability with the true internal variability obtained from the debug output. The RMSE provides a quantitative measure of the agreement between the two time series.

Run the script:

python fig-timeseries_internal_variability.py

Then repeat the analysis with a larger number of years:

N_YEAR = 15

Questions:

  • How well does the estimated internal variability match the true internal variability?
  • Does the RMSE decrease when N_YEAR increases?
  • Why does removing the climatological seasonal cycle help isolate internal variability?
  • Why is the estimated internal variability not exactly identical to the true internal variability?

Key Points

  • The Lorenz model provides a simple framework for investigating internal climate variability.
  • A periodic external forcing can represent the seasonal cycle in an idealized climate model.
  • The total climate state can be interpreted as the sum of the externally forced response and internally generated variability.
  • The climatological seasonal cycle is estimated by averaging many years of model output.
  • Averaging over a larger number of years suppresses internal variability and improves the estimate of the externally forced seasonal cycle.
  • Internal variability can be estimated by removing the climatological seasonal cycle from the total climate state.
  • Comparing the estimated components with the debug output provides a quantitative assessment of the decomposition using RMSE.

Homework

Repeat the seasonal-cycle experiment using a longer integration.

Modify the model integration period in run_model_climate-v3.sh:

TEND="420"

which corresponds to a 30-year simulation following the spin-up.

Run the model again:

./run_model_climate-v3.sh

Then update the plotting scripts to use the longer simulation, for example,

TSTR = 60
TEND = 420

N_YEAR = 30

Generate the following two figures:

  1. Estimated externally forced seasonal cycle (fig-seasonal_cycle-rmse.py)
  • Plot the individual years together with the climatological seasonal cycle.
  • Compare the estimated climatology with the true externally forced component.
  • Display the RMSE and N_YEAR in each panel.
  1. Estimated internal variability (fig-timeseries_internal_variability.py)
  • Plot the estimated internal variability together with the true internal variability.
  • Display the RMSE and N_YEAR in each panel.

Submit:

  • fig-seasonal_cycle-rmse.png
  • fig-timeseries_internal_variability.png