7.3 External Forcing II: Climate Model Projection

Learning Objectives

After completing this lecture, students should be able to:

explain the difference between climate variability and climate change
perform pre-industrial control and forced climate model experiments
understand why multiple ensemble simulations are required for climate projections
estimate externally forced climate change using the ensemble mean
estimate internal variability from deviations about the ensemble mean

Concept

The previous section introduced a periodically varying external forcing to represent the seasonal cycle. Although the forcing repeated every model year, the climate state varied from year to year because of internally generated variability.

Climate change can be represented in a similar framework by prescribing a persistent external forcing rather than a periodic one. For example, increasing greenhouse gas concentrations produce a long-term radiative forcing that shifts the climate toward a new equilibrium.

In this section, we prescribe a constant external forcing. The resulting climate response consists of two components:

a long-term forced response caused by the external forcing,
internally generated variability arising from the chaotic dynamics of the Lorenz model.

Because internal variability differs among simulations, a single model realization cannot perfectly represent the externally forced climate change. Instead, climate projections are commonly performed using an ensemble of simulations that begin from slightly different initial conditions.

Averaging across the ensemble suppresses the internally generated variability and provides an estimate of the externally forced climate response. The remaining deviations from the ensemble mean represent internal variability.

This decomposition forms the basis of modern climate projection studies using Earth system models.

Model / Equation

As in the previous sections, the total climate state is written as

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

The internally generated component evolves according to the Lorenz equations,

\[ \frac{dX_{\rm int}}{dt} = F_{\rm int}(X_{\rm int}), \]

while the externally forced component is prescribed by a constant forcing,

\[ \frac{dX_{\rm ext}}{dt} = F_{\rm ext}, \]

where

\[ F_{\rm ext} = (F_x,F_y,F_z) \]

is constant throughout the simulation.

Unlike the seasonal-cycle experiment, the forcing does not vary with time. Consequently, the externally forced component gradually shifts the mean climate state, producing a climate trend analogous to long-term global warming.

To estimate the externally forced response, an ensemble of simulations is performed. The ensemble mean is

\[ \overline{X} = \frac{1}{N} \sum_{i=1}^{N} X_i, \]

where \(N\) is the number of ensemble members.

The internal variability for each ensemble member is estimated as

\[ X_i' = X_i - \overline{X}. \]

As the number of ensemble members increases, the ensemble mean provides an increasingly accurate estimate of the externally forced climate response.

Sample Code

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

Model_src-v3.0.tar.gz

Download the archive and extract it using

tar -xzvf Model_src-v3.1.tar.gz

This creates the directory

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

The individual source files are also available below.

run_model_climate-v3.sh — Configure the climate model experiment.
Lmodel_climate-v3.py — Lorenz climate model with constant external forcing.
fig-timeseries_climate.py — Plot the total climate state.
fig-timeseries_climate_debug.py — Plot the total, internal, and externally forced components for debugging.

Compared with the previous lecture, only the external forcing has been modified from a periodic forcing to a constant forcing. The model structure, numerical integration, and visualization scripts remain unchanged. This allows us to focus on how different types of external forcing produce different climate responses while using the same Lorenz climate model framework.

Exercise

Step 1. Model Setup and Spin-Up

Create a new working directory for this exercise:

mkdir Model_EXP7_3
cd Model_EXP7_3

Copy the source code:

cp ../Model_src-v3/*.* .
mkdir Spinup

Move to the spin-up directory and copy the source code:

cd Spinup
cp ../*.* ./

Verify that the directory contains

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

Edit run_model_climate-v3.sh and prescribe

REST="False"
TSTR="0"
TEND="100"

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

DEBUG_OUTPUT="False"

Run the model:

sh run_model_climate-v3.sh

The restart file generated at the end of the integration will be used to initialize the following experiments.

Step 2. Pre-industrial Control Simulation (piControl)

Return to the experiment directory, then set up the environment to restart from the spin-up:

cd ..

mkdir piControl
cd piControl

cp ../*.* ./

mkdir rest

Copy the restart file from the spin-up experiment:

cp ../Spinup/rest/restart.0100.nc ./rest/

Edit run_model_climate-v3.sh:

REST="True"

TSTR="100"
TEND="200"

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

DEBUG_OUTPUT="False"

Run the model:

sh ./run_model_climate-v3.sh

Visualize the result:

mkdir figures

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

cd figures

Open fig-timeseries_climate.py and set

TSTR = 100
TEND = 200

Then run

python fig-timeseries_climate.py

Questions

Does the climate remain statistically stationary?
Is the variability entirely internally generated?

Step 3. Climate Projection (ENS001)

Return to the experiment directory and create the first ensemble member:

cd ..

mkdir ENS001
cd ENS001

cp ../*.* ./

mkdir rest

Copy the restart file from the end of the pre-industrial control simulation:

cp ../piControl/rest/restart.0110.nc  ./rest/restart.0110.piControl.nc

This experiment will begin at model year 300. First, convert the restart file into the CDL (network Common Data form Language) text format:

ncdump rest/restart.0110.piControl.nc > ./restart.0300.cdl

Open restart.0300.cdl using a text editor and locate the line

time = 110 ;

Change it to

time = 300 ;

Then regenerate the NetCDF restart file:

ncgen -o rest/restart.0300.nc  ./restart.0300.cdl

Verify that the restart file has been updated correctly:

ncdump rest/restart.0300.nc

Edit run_model_climate-v3.sh:

REST="True"

TSTR="300"
TEND="400"

AMP_X="0.05"
AMP_Y="0.1"
AMP_Z="0.15"

DEBUG_OUTPUT="True"

Run the model:

sh run_model_climate-v3.sh

Visualize the result:

mkdir figures

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

cd figures

Open fig-timeseries_climate.py and set

TSTR = 300
TEND = 400

Run the plotting script:

python fig-timeseries_climate.py

Questions

Does the model exhibit a long-term climate trend under constant external forcing?
Is internal variability still present?
Why does a single ensemble member not perfectly represent the externally forced climate response?

Step 4. Additional Ensemble Members

Repeat Step 3 to generate additional ensemble members. Each member uses the same model configuration and the same external forcing, but begins from a different initial condition obtained from the pre-industrial control simulation.

Create the second ensemble member by copying the first experiment:

cd ..
mkdir ENS002

cp ENS001/Lmodel_climate-v3.py ENS001/run_model_climate-v3.sh ENS002/

cd ENS002

Replace the restart file with a different initial condition from the pre-industrial control simulation:

mkdir rest
cp ../piControl/rest/restart.0120.nc ./rest/restart.0120.piControl.nc

Convert the restart file to CDL format:

ncdump rest/restart.0120.piControl.nc > restart.0300.cdl

Open restart.0300.cdl and change

time = 120 ;

to

time = 300 ;

Then regenerate the restart file:

ncgen -o rest/restart.0300.nc restart.0300.cdl

Verify the restart file:

ncdump rest/restart.0300.nc

Since the entire experiment was copied from ENS001, no further modifications to run_model_climate-v3.sh are required. The external forcing and all model parameters remain identical.

Run the model:

sh run_model_climate-v3.sh

Repeat this procedure to create additional ensemble members. For example,

ENS003  ← restart.0130.nc
ENS004  ← restart.0140.nc
ENS005  ← restart.0150.nc

Only the restart file should differ among the ensemble members.

Visualize each ensemble member using the same plotting script from Step 3.

Questions

Which features are common among all ensemble members?
Which features differ from member to member?
Why do the simulations diverge even though the external forcing is identical?
Why is it important that all ensemble members use exactly the same model configuration except for the initial condition?

Step 5. Visualize Ensemble Members

Return to the main experiment directory:

cd ..

Create a directory for ensemble visualization:

mkdir figures

Copy the time-series plotting script and rename it:

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

Move to the figures directory:

cd figures

Open fig-timeseries_ensemble.py and replace the single-member file-reading section with the following code.

First, begin with two ensemble members:

# ============================================================
# Figure settings
# ============================================================
script_name = os.path.splitext(os.path.basename(__file__))[0]
outfile_png = script_name + ".png"

ENS = [1, 2]

ENSEMBLE_MEMBERS = [
    f"ENS{member:03d}" for member in ENS
]

TSTR = 300
TEND = 400

Then read all selected ensemble members:

# ============================================================
# Read model output
# ============================================================
Xens = []
Yens = []
Zens = []

for member in ENSEMBLE_MEMBERS:

    fin = xr.open_mfdataset(
        f"../{member}/out/out.*.nc",
        combine="by_coords"
    )

    Xens.append(fin.X.sel(time=slice(TSTR, TEND)).values)
    Yens.append(fin.Y.sel(time=slice(TSTR, TEND)).values)
    Zens.append(fin.Z.sel(time=slice(TSTR, TEND)).values)

    if member == ENSEMBLE_MEMBERS[0]:
        time = fin.time.sel(time=slice(TSTR, TEND)).values

    fin.close()

Xens = np.asarray(Xens)
Yens = np.asarray(Yens)
Zens = np.asarray(Zens)

Replace the original plotting commands such as

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

with

for iens, member in enumerate(ENSEMBLE_MEMBERS):

    ax[0].plot( time, Xens[iens, :], linewidth=1.0, label=member)
    ax[1].plot( time, Yens[iens, :], linewidth=1.0, label=member)
    ax[2].plot( time, Zens[iens, :], linewidth=1.0, label=member)

Add a legend:

ax[0].legend(loc="best", fontsize=8)

Run the script:

python fig-timeseries_ensemble.py

After confirming that the two-member figure works, increase the number of ensemble members:

ENS = [1, 2, 3, 4, 5]

Run the script again:

python fig-timeseries_ensemble.py

Questions:

Do the ensemble members show exactly the same trajectory?
Which features are common among the ensemble members?
Which features differ among ensemble members?
Why do the ensemble members differ even though the same external forcing is prescribed?

Step 6. Estimate the Ensemble Mean

The ensemble mean provides an estimate of the externally forced climate response because internally generated variability differs among ensemble members and is reduced by averaging.

Create a new plotting script:

cp fig-timeseries_ensemble.py fig-timeseries_ensMEAN.py

Open fig-timeseries_ensMEAN.py.

First, use the same ensemble-member reading section from Step 5:

ENS = [1, 2]

Calculate the ensemble mean after read the data:

Xmean = np.mean(Xens, axis=0)
Ymean = np.mean(Yens, axis=0)
Zmean = np.mean(Zens, axis=0)

Plot each ensemble member using thin, light-gray lines:

for iens in range(len(ENSEMBLE_MEMBERS)):

    ax[0].plot(time, Xens[iens, :], color="lightgray", linewidth=0.8)
    ax[1].plot(time, Yens[iens, :], color="lightgray", linewidth=0.8)
    ax[2].plot(time, Zens[iens, :], color="lightgray", linewidth=0.8)

Overlay the ensemble mean using a thick colored line:

ax[0].plot(time, Xmean, color="C0", linewidth=3.0, label="Ensemble Mean")
ax[1].plot(time, Ymean, color="C0", linewidth=3.0, label="Ensemble Mean")
ax[2].plot(time, Zmean, color="C0", linewidth=3.0, label="Ensemble Mean")

Run the script:

python fig-timeseries_ensMEAN.py

After confirming that the ensemble mean is plotted correctly, add the true externally forced component from the debug output.

Inside the ensemble-member reading loop, add:

Xext_ens = []
Yext_ens = []
Zext_ens = []

and read the externally forced components:

Xext_ens.append(fin.Xext.sel(time=slice(TSTR, TEND)).values)
Yext_ens.append(fin.Yext.sel(time=slice(TSTR, TEND)).values)
Zext_ens.append(fin.Zext.sel(time=slice(TSTR, TEND)).values)

After the loop, convert the externally forced components to NumPy arrays:

Xext_ens = np.asarray(Xext_ens)
Yext_ens = np.asarray(Yext_ens)
Zext_ens = np.asarray(Zext_ens)

Because the externally forced component is identical for all ensemble members, use the first ensemble member as the reference:

Xext = Xext_ens[0, :]
Yext = Yext_ens[0, :]
Zext = Zext_ens[0, :]

Adjust the mean of the externally forced component so that it can be compared directly with the ensemble mean:

Xext = Xext + (np.mean(Xmean) - np.mean(Xext))
Yext = Yext + (np.mean(Ymean) - np.mean(Yext))
Zext = Zext + (np.mean(Zmean) - np.mean(Zext))

Calculate the root-mean-square error (RMSE) between the ensemble mean and the true externally forced component:

rmseX = np.sqrt(np.mean((Xmean - Xext)**2))
rmseY = np.sqrt(np.mean((Ymean - Yext)**2))
rmseZ = np.sqrt(np.mean((Zmean - Zext)**2))

Overlay the true externally forced component using a dashed black line:

ax[0].plot(time, Xext, color="k", linestyle="--", linewidth=2.0,
           label="True External")
ax[1].plot(time, Yext, color="k", linestyle="--", linewidth=2.0,
           label="True External")
ax[2].plot(time, Zext, color="k", linestyle="--", linewidth=2.0,
           label="True External")

Display the ensemble size and RMSE in each panel:

ax[0].text(
    0.02, 0.88,
    f"N_ENS = {len(ENSEMBLE_MEMBERS)}\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_ENS = {len(ENSEMBLE_MEMBERS)}\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_ENS = {len(ENSEMBLE_MEMBERS)}\nRMSE = {rmseZ:.2f}",
    transform=ax[2].transAxes,
    bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)

Run the script again:

python fig-timeseries_ensMEAN.py

Increase the number of ensemble members and repeat the analysis:

ENS = [1, 2, 3, 4, 5]

Questions

How does the ensemble mean change as the number of ensemble members increases?
How does the RMSE change as additional ensemble members are included?
Why is the ensemble mean smoother than the individual ensemble members?
Why does the ensemble mean provide a better estimate of the externally forced climate response?

Key Points

A spin-up simulation is used to obtain an equilibrated initial climate state before conducting scientific experiments.
A pre-industrial control (piControl) simulation represents natural climate variability in the absence of external forcing.
A climate projection is performed by applying external forcing to the equilibrated climate state.
Ensemble members are generated by using different initial conditions while keeping the model configuration and external forcing identical.
The spread among ensemble members represents internal climate variability.
The ensemble mean suppresses internal variability and provides an estimate of the externally forced climate response.
Increasing the number of ensemble members generally reduces the error in the estimated externally forced component.

Homework

Repeat the climate projection experiment using 10 ensemble members.

Use restart files from the pre-industrial control simulation, for example,

restart.0110.nc
restart.0120.nc
...
restart.0200.nc

to create

ENS001
ENS002
...
ENS010

All ensemble members should use the same external forcing and the same model configuration. Only the initial condition should differ.

Modify fig-timeseries_ensMEAN.py to use

ENS = [1,2,3,4,5,6,7,8,9,10]

Generate the final figure showing

individual ensemble members as thin light-gray lines,
the ensemble mean as a thick colored line,
the true externally forced component as a dashed black line,
the RMSE between the ensemble mean and the true externally forced component.

Submit:

The final figure comparing the 10-member ensemble mean with the true externally forced component.