Chapter 9.2 Ensemble Forecasting II: Ensemble Size and Predictability

Learning Objectives

After completing this chapter, you should be able to

Develop an automated framework to perform ensemble hindcast experiments at multiple initialization times.
Organize ensemble forecasts into multidimensional arrays for forecast verification.
Calculate the ensemble mean from multiple ensemble members.
Reuse deterministic forecast-verification code to evaluate ensemble forecasts.
Evaluate forecast skill using the Anomaly Correlation Coefficient (ACC), Root Mean Square Error (RMSE), and RMSE Skill Score (RMSESS).
Investigate how increasing ensemble size influences forecast skill.
Discuss the relationship between ensemble size, forecast uncertainty, and predictability.

Concepts

In Chapter 9.1, ensemble forecasting was introduced by generating perturbed initial conditions and visualizing the ensemble mean and spread for a single forecast. In this chapter, the same framework is extended to evaluate predictive skill using multiple hindcast experiments.

Rather than evaluating a single initialization, ensemble forecasts are generated for several initialization times. Each hindcast experiment contains multiple ensemble members, allowing the ensemble mean to be compared directly with the deterministic forecast.

The ensemble forecasts are organized into a three-dimensional array,

Xens(ensemble, experiment, lead)
Yens(ensemble, experiment, lead)
Zens(ensemble, experiment, lead)

where

ensemble represents the ensemble member,
experiment represents the hindcast initialization time, and
lead represents the forecast lead time.

The ensemble mean is then calculated by averaging over the ensemble dimension,

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

Because these arrays have the same dimensions as those used in Chapter 8.3, the existing forecast-verification framework can be reused without modification.

This chapter also investigates one of the most important questions in ensemble prediction:

How many ensemble members are needed to improve forecast skill?

Increasing the ensemble size generally produces a more stable ensemble mean by reducing the influence of individual perturbations. However, larger ensembles also require greater computational resources, and the improvement in predictive skill often becomes smaller as additional members are included.

Model and Equations

Suppose that an ensemble forecast consists of \(N_{\mathrm{ens}}\) ensemble members. For each hindcast experiment and forecast lead time, the forecast from the \(i\)-th ensemble member is denoted by

\[ X_i(e,l), \]

where

\(e\) denotes the hindcast experiment (initialization time),
\(l\) denotes the forecast lead time, and
\(i=1,2,\ldots,N_{\mathrm{ens}}\).

Similarly, ensemble forecasts are generated for \(Y\) and \(Z\).

Ensemble Mean

The ensemble mean is computed by averaging over all ensemble members,

\[ X_f(e,l) = \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} X_i(e,l). \]

The same calculation is applied to the \(Y\) and \(Z\) components,

\[ Y_f(e,l) = \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} Y_i(e,l), \]

\[ Z_f(e,l) = \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} Z_i(e,l). \]

The ensemble forecasts are organized as

Xens(ensemble, experiment, lead)
Yens(ensemble, experiment, lead)
Zens(ensemble, experiment, lead)

After averaging over the ensemble dimension, the resulting arrays become

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

These arrays have the same dimensions as the deterministic hindcast forecasts developed in Chapter 8.3 and therefore can be verified using the same forecast-skill calculations.

Anomaly Correlation Coefficient (ACC)

The Anomaly Correlation Coefficient is defined as

\[ \mathrm{ACC} = \frac{ \sum_{e} \left(X_f-\overline{X_f}\right) \left(X_{\mathrm{true}}-\overline{X_{\mathrm{true}}}\right) }{ \sqrt{ \sum_{e} \left(X_f-\overline{X_f}\right)^2 } \sqrt{ \sum_{e} \left(X_{\mathrm{true}}-\overline{X_{\mathrm{true}}}\right)^2 } }. \]

ACC measures the linear correspondence between the forecast and the true state.

Root Mean Square Error (RMSE)

The Root Mean Square Error is

\[ \mathrm{RMSE} = \sqrt{ \frac{1}{N_{\mathrm{exp}}} \sum_{e=1}^{N_{\mathrm{exp}}} \left( X_f-X_{\mathrm{true}} \right)^2 }, \]

where \(N_{\mathrm{exp}}\) is the number of hindcast experiments.

RMSE measures the magnitude of the forecast error.

RMSE Skill Score (RMSESS)

The RMSE Skill Score compares the forecast with a climatological reference forecast,

\[ \mathrm{RMSESS} = 1- \frac{\mathrm{RMSE}} {\mathrm{RMSE}_{\mathrm{climatology}}}. \]

Its interpretation is

RMSESS	Interpretation
1	Perfect forecast
\(>0\)	Better than climatology
0	Equal to climatology
\(<0\)	Worse than climatology

As the ensemble size increases, the ensemble mean generally becomes more stable because random perturbations tend to cancel one another. Consequently, ensemble forecasts often produce larger ACC, smaller RMSE, and larger RMSE Skill Score than a single deterministic forecast.

Sample Code

The sample programs in this section extend the ensemble-forecast framework developed in Chapter 9.1. Multiple hindcast experiments are generated automatically, and the resulting ensemble forecasts are organized into multidimensional arrays for forecast verification.

Download the source code archive,

Model_src-v3.0.tar.gz

and extract it using

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

This creates the directory

Model_src-v3.0/
├── figures/
│   ├── fig-timeseries_hindcast_ENS_acc-template.py
│   └── fig-hindcast_skill-template.py
├── utils/
│   ├── get_observation.py
│   ├── create_initial_obs.py
│   ├── create_initial_obs-perturbed.py
│   └── run_hindcast.sh
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

The individual source files used in this chapter are listed below.

File	Description
`run_model_climate-v3.sh`	Defines the experiment settings and executes the Lorenz climate model.
`Lmodel_climate-v3.py`	Performs numerical integration of the Lorenz climate model.
`create_initial_obs-perturbed.py`	Creates deterministic and perturbed initial conditions from synthetic observations.
`run_hindcast.sh`	Automatically generates the deterministic hindcast experiments (C00) for multiple initialization times.
`fig-timeseries_hindcast_ENS_acc-template.py`	Visualizes ensemble forecasts and compares the true state, individual ensemble members, and the ensemble mean.
`fig-hindcast_skill-template.py`	Calculates ACC, RMSE, and RMSE Skill Score for deterministic and ensemble-mean forecasts.

The hindcast experiments are initialized at

T_INITIAL = 55, 60, 65, 70, 75, 80, 85, and 90

For each initialization time, the deterministic forecast (C00) and multiple perturbed ensemble members are stored in separate directories.

For example,

HINDCAST/
├── 0055/
│   ├── C00/
│   ├── P01/
│   ├── M01/
│   ├── P02/
│   ├── M02/
│   ├── P03/
│   └── M03/
├── 0060/
├── 0065/
├── ...
└── run_ENS-v1.sh

The ensemble forecasts are organized into three-dimensional arrays,

Xens(ensemble, experiment, lead)
Yens(ensemble, experiment, lead)
Zens(ensemble, experiment, lead)

The ensemble mean is then calculated as

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

which have the same dimensions as the deterministic forecasts developed in Chapter 8.3.

Consequently, the existing forecast-verification framework can be reused directly to calculate

Anomaly Correlation Coefficient (ACC),
Root Mean Square Error (RMSE), and
RMSE Skill Score (RMSESS).

By modifying only

ENSEMBLE_MEMBERS

the same verification program can evaluate different ensemble sizes, automatically updating the figure title and output filename. This modular design allows the influence of ensemble size on forecast skill to be investigated with minimal changes to the source code.

Exercise

Step 1. Generate the Control Simulation and Synthetic Observations

In this exercise, we first generate the control (CTL) simulation, which represents the true evolution of the Lorenz climate system. The control simulation serves as the reference ("truth") for evaluating the deterministic and ensemble hindcasts developed later in this chapter.

Create the working directories,

mkdir Chapter9_2
cd Chapter9_2

mkdir CTL
cd CTL

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

The directory should initially contain

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

Run the Control Simulation

Open run_model_climate-v3.sh and configure the experiment as

REST=False
TSTR=0
TEND=100

Then execute

bash run_model_climate-v3.sh

The model automatically creates the required output directories and integrates the control simulation over the period \(t=0\)–100.

Generate Synthetic Observations

The deterministic and ensemble hindcasts will be initialized from synthetic observations. In this chapter, use a smaller observation error than in Chapter 9.1,

OERR = 1.0

to represent a more accurate observing system.

Copy the observation-generation script,

cp ../../Model_src-v3.0/utils/get_observation.py .

Modify the following settings,

OERR = 1.0

TSTR = 50
TEND = 100

OBS_TSTR = 55.0
OBS_TEND = 95.0

and execute

python get_observation.py

This creates

observation.nc

which will be used to initialize all hindcast experiments in this chapter.

Visualize and Verify the CTL Simulation

Return to the directory above and create a directory for the figure scripts,

cd ../

mkdir Figure
cd Figure

Copy the observation figure template,

cp ../../../Model_src-v3.0/figures/fig-timeseries_obs.py .

Open

fig-timeseries_obs.py

and modify

FDIR = "../CTL/out/"
OBSFILE = "../CTL/observation.nc"

TSTR = 50
TEND = 100

Run the figure script,

python fig-timeseries_obs.py

Verify that the synthetic observations fluctuate around the true control simulation because of the prescribed observation error.

Verify the Results

After completing this step, the directory should resemble

CTL/
├── out/
│   ├── out.0000.nc
│   ├── out.0001.nc
│   └── ...
├── rest/
├── observation.nc
├── get_observation.py
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

Figure/
└── fig-timeseries_obs.py

The CTL simulation represents the true climate evolution, while observation.nc provides the synthetic observations that will be used to initialize the deterministic and ensemble hindcast experiments developed in the following exercises.

Step 2. Automate the Deterministic Hindcast Experiments (C00)

In this step, we generate a series of deterministic hindcast experiments initialized directly from the synthetic observations. These experiments are referred to as C00 because the perturbation amplitude is set to zero.

The hindcasts are initialized every five time units,

T_INITIAL = 55, 60, 65, 70, 75, 80, 85, and 90

and each hindcast is integrated for five time units.

Create the Hindcast Directory

Return to the Chapter9_2 directory,

cd ../

Create the working directory,

mkdir HINDCAST
cd HINDCAST

Copy the perturbed initial-condition generator,

cp ../../Model_src-v3.0/utils/create_initial_obs-perturbed.py .

Before generating the restart file, update the observation file path in

create_initial_obs-perturbed.py

Because the script will be executed under

HINDCAST/TTTT/C00/

the relative path to the observation file must be changed to

OBSFILE = "../../../CTL/observation.nc"

This modification ensures that the initialization program reads the synthetic observations generated in the CTL experiment.

The directory should initially contain

HINDCAST/
└── create_initial_obs-perturbed.py

Create the Automated Hindcast Script

Create a shell script,

vi run_C00.sh

and add the following code.

#!/bin/bash

# ============================================================
# Automated Deterministic Hindcast Experiments (C00)
# ============================================================

# Stop the script immediately if any command fails.
set -e

# ============================================================
# Experiment settings
# ============================================================
T_INITIAL_LIST=(55 60 65 70 75 80 85 90)

FORECAST_LENGTH=5

SEED=1
AMP=0.0


# ============================================================
# Run hindcast experiments
# ============================================================
for T_INITIAL in "${T_INITIAL_LIST[@]}"
do

    T_END=$((T_INITIAL + FORECAST_LENGTH))

    EXP_NAME=$(printf "%04d" "${T_INITIAL}")
    RUN_SCRIPT="run_model_climate-v3_${EXP_NAME}.sh"

    echo "========================================"
    echo "Create C00 Hindcast Experiment"
    echo "========================================"
    echo "T_INITIAL : ${T_INITIAL}"
    echo "T_END     : ${T_END}"
    echo "Directory : ${EXP_NAME}/C00"
    echo "========================================"

    mkdir -p "${EXP_NAME}/C00"
    cd "${EXP_NAME}/C00"

    # Copy the model source code
    cp ../../../CTL/run_model_climate-v3.sh .
    cp ../../../CTL/Lmodel_climate-v3.py .

    # Copy the initialization program
    cp ../../create_initial_obs-perturbed.py .

    # Create the restart file
    python create_initial_obs-perturbed.py \
        --t_initial="${T_INITIAL}" \
        --seed="${SEED}" \
        --amp="${AMP}"

    # Configure the hindcast experiment
    sed \
        -e 's/^REST=.*/REST="True"/' \
        -e "s/^TSTR=.*/TSTR=\"${T_INITIAL}\"/" \
        -e "s/^TEND=.*/TEND=\"${T_END}\"/" \
        run_model_climate-v3.sh \
        > "${RUN_SCRIPT}"

    echo "----------------------------------------"
    echo "Model Settings"
    echo "----------------------------------------"
    grep -E '^(REST|TSTR|TEND)=' "${RUN_SCRIPT}"
    echo "----------------------------------------"

    # Run the model
    bash "${RUN_SCRIPT}"

    cd ../../

    echo "========================================"
    echo "C00 experiment ${EXP_NAME} completed."
    echo "========================================"

done

echo "========================================"
echo "All C00 hindcast experiments completed."
echo "========================================"

This script automatically performs the following tasks for each initialization time:

Creates a hindcast directory for the corresponding initialization time.
Creates a C00 experiment directory.
Copies the model source code from the CTL experiment.
Creates the restart file using the synthetic observation with

AMP = 0.0

so that no perturbation is added to the observed initial condition.

Configures the deterministic hindcast experiment.
Executes the model simulation.

Run the C00 Hindcast Experiments

Make the script executable,

chmod +x run_C00.sh

Run all deterministic hindcasts,

./run_C00.sh

The script generates the following directory structure,

HINDCAST/
├── 0055/
│   └── C00/
├── 0060/
│   └── C00/
├── 0065/
│   └── C00/
├── 0070/
│   └── C00/
├── 0075/
│   └── C00/
├── 0080/
│   └── C00/
├── 0085/
│   └── C00/
├── 0090/
│   └── C00/
├── create_initial_obs-perturbed.py
└── run_C00.sh

Each experiment directory contains a deterministic hindcast initialized directly from the synthetic observation.

For example,

HINDCAST/
└── 0060/
    └── C00/
        ├── out/
        ├── rest/
        │   └── restart.0060.nc
        ├── create_initial_obs-perturbed.py
        ├── Lmodel_climate-v3.py
        ├── run_model_climate-v3.sh
        └── run_model_climate-v3_0060.sh

Visualize and Verify the C00 Hindcasts

Return to the Figure directory,

cd ../../Figure

Copy the hindcast-verification script,

cp ../../../Model_src-v3.0/figures/fig-timeseries_hindcast_acc.py \
   fig-timeseries_C00_acc.py

Open

fig-timeseries_C00_acc.py

and revise the figure description and title so that they refer to the C00 hindcasts.

For example,

figure_title = "Lorenz Model: ACC of X for C00 Hindcasts"

Update the hindcast output directory to reflect the directory structure developed in this chapter,

FDIR_CTL = "../CTL/out/"
FDIR_HINDCAST = "../HINDCAST/dddd/C00/out/"

Continue to use the initialization times

HINDCAST_EXPERIMENTS = [
    "0055",
    "0060",
    "0065",
    "0070",
    "0075",
    "0080",
    "0085",
    "0090",
]

and the forecast lead times

LEAD_TIMES = [0, 1, 2, 3, 4]

The following code replaces dddd with each four-digit initialization time,

FDIR_PRD = FDIR_HINDCAST.replace( "dddd", experiment)

For example,

experiment = "0060"

produces

../HINDCAST/0060/C00/out/

Revise the hindcast labels so that the forecasts are identified as C00 hindcasts. For example,

if iexp == 0 and ilead == 0:
    hindcast_label = "C00 hindcasts"
else:
    hindcast_label = None

You may also change

forecast_label = "Forecast state"

to

forecast_label = "C00 state"

Run the figure script,

python fig-timeseries_C00_acc.py

The figure displays five panels corresponding to forecast lead times

0, 1, 2, 3, and 4

Each panel shows

the CTL simulation, representing the true evolution,
the eight C00 hindcast trajectories,
the true state at the selected verification time,
the corresponding C00 forecast state, and
the ACC calculated across all initialization experiments.

Simplify the C00 Hindcast Figure

To prepare for the ensemble experiments in the next step, simplify the figure by displaying only the verification states rather than the complete forecast trajectories.

Comment out the commands that plot the C00 forecast trajectory,

# ax[ilead].plot(
#     time_prd,
#     Xprd,
#     color="C1",
#     linewidth=1.0,
#     alpha=0.8,
#     label=hindcast_label
# )

Also comment out the commands that connect the true and forecast states,

# ax[ilead].plot(
#     [verification_time, verification_time],
#     [true_value, forecast_value],
#     color="0.5",
#     linewidth=0.8
# )

Retain only

the CTL trajectory,
the true state (blue circle), and
the C00 forecast state (red ×).

Run the revised figure script,

python fig-timeseries_C00_acc.py

The resulting figure provides a clean visualization of the deterministic hindcast verification and serves as the foundation for the ensemble verification figure developed in the next exercise.

Verify the Results

Confirm that the script reads the following experiments:

HINDCAST/0055/C00/
HINDCAST/0060/C00/
HINDCAST/0065/C00/
HINDCAST/0070/C00/
HINDCAST/0075/C00/
HINDCAST/0080/C00/
HINDCAST/0085/C00/
HINDCAST/0090/C00/

Examine how ACC changes with forecast lead time.

Consider the following questions.

Is ACC highest at short forecast lead times?
Does ACC decrease as forecast lead time increases?
Are all C00 hindcasts equally accurate?
How does observational error in the initial condition influence the C00 forecasts?

The C00 skill provides the deterministic baseline that will later be compared with ensemble-mean forecasts of different ensemble sizes.

Step 3. Generate Three-Member Ensemble Hindcasts

In this step, we extend the deterministic C00 hindcasts by generating one symmetric perturbation pair, P01 and M01, for each initialization time.

Together,

C00
P01
M01

form a three-member ensemble hindcast.

The positively and negatively perturbed members use the same random seed and opposite perturbation amplitudes,

P01    SEED = 1    AMP = +1.0
M01    SEED = 1    AMP = -1.0

This symmetric construction ensures that the ensemble mean is equal to the observed initial condition at the forecast initialization time.

Create an Automated Ensemble-Hindcast Script

Return to the HINDCAST directory,

cd ../HINDCAST

Copy the deterministic hindcast script,

cp run_C00.sh run_ENS-v1.sh

Open

run_ENS-v1.sh

and replace the experiment settings with

# ============================================================
# Experiment settings
# ============================================================
T_INITIAL_LIST=(55 60 65 70 75 80 85 90)

MEMBERS=("P01" "M01")
SEEDS=(1 1)
AMPS=(1.0 -1.0)

Replace the loop over the C00 experiments with the following loop.

#!/bin/bash

# ============================================================
# Automated Three-Member Ensemble Hindcast Experiments
#
# Existing member:
#   C00
#
# Members generated by this script:
#   P01
#   M01
# ============================================================

# Exit immediately if a command fails.
set -e


# ============================================================
# Experiment settings
# ============================================================
T_INITIAL_LIST=(55 60 65 70 75 80 85 90)

MEMBERS=("P01" "M01")
SEEDS=(1 1)
AMPS=(1.0 -1.0)


# ============================================================
# Generate ensemble hindcasts
# ============================================================
for T_INITIAL in "${T_INITIAL_LIST[@]}"
do

    EXP_NAME=$(printf "%04d" "${T_INITIAL}")

    for i in "${!MEMBERS[@]}"
    do

        MEMBER=${MEMBERS[$i]}
        SEED=${SEEDS[$i]}
        AMP=${AMPS[$i]}

        RUN_SCRIPT="run_model_climate-v3_${EXP_NAME}.sh"

        echo "========================================"
        echo "Create Ensemble Hindcast"
        echo "========================================"
        echo "Initialization : ${EXP_NAME}"
        echo "Member         : ${MEMBER}"
        echo "Seed           : ${SEED}"
        echo "Amplitude      : ${AMP}"
        echo "========================================"

        mkdir -p "${EXP_NAME}/${MEMBER}"
        cd "${EXP_NAME}/${MEMBER}"

        # Copy the verified model configuration
        cp ../C00/${RUN_SCRIPT} .
        cp ../C00/Lmodel_climate-v3.py .

        # Copy the perturbed initial-condition generator
        cp ../../create_initial_obs-perturbed.py .

        # Generate the perturbed restart file
        python create_initial_obs-perturbed.py \
            --t_initial="${T_INITIAL}" \
            --seed="${SEED}" \
            --amp="${AMP}"

        # Run the hindcast
        bash ${RUN_SCRIPT}

        cd ../../

        echo "${EXP_NAME}/${MEMBER} completed."

    done

done

echo "========================================"
echo "Three-member ensemble hindcasts completed."
echo "========================================"

Because each ensemble member copies the model source code from the corresponding C00 experiment, all members use exactly the same model configuration. The only difference among the ensemble members is the perturbed initial condition.

Run the Ensemble Hindcasts

Make the script executable,

chmod +x run_ENS-v1.sh

Run the ensemble experiments,

./run_ENS-v1.sh

The resulting directory structure should resemble

HINDCAST/
├── 0055/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0060/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0065/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0070/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0075/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0080/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0085/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── 0090/
│   ├── C00/
│   ├── P01/
│   └── M01/
├── create_initial_obs-perturbed.py
├── run_C00.sh
└── run_ENS-v1.sh

Verify the Results

Confirm that each initialization experiment contains

C00
P01
M01

For example,

HINDCAST/
└── 0060/
    ├── C00/
    │   ├── out/
    │   └── rest/restart.0060.nc
    ├── P01/
    │   ├── out/
    │   └── rest/restart.0060.nc
    └── M01/
        ├── out/
        └── rest/restart.0060.nc

The positive and negative members use the same random seed but opposite perturbation amplitudes,

P01 : seed = 1, amp = +1.0
M01 : seed = 1, amp = -1.0

Therefore, the perturbations are symmetric about the observed initial condition, and the three-member ensemble provides the first estimate of forecast uncertainty.

In the next step, the deterministic verification figure will be extended to calculate the ensemble mean, ensemble spread, and predictive skill for the three-member ensemble.

Step 4. Visualize the Three-Member Ensemble Hindcasts

In this step, we visualize the three-member ensemble hindcasts generated in the previous exercise. We first use the provided template to display the individual ensemble members,

C00
P01
M01

at each verification time. We then modify the same figure to summarize the forecasts using the ensemble mean and ensemble spread.

The forecast data are organized as

Xens(ensemble, initial, lead)

where

ensemble identifies the ensemble member,
initial identifies the forecast initialization time, and
lead identifies the forecast lead time.

Visualize the Individual Ensemble Members

Return to the Figure directory,

cd ../Figure

Copy the ensemble-verification template,

cp ../../../Model_src-v3.0/figures/fig-timeseries_hindcast_ENS_acc-template.py \
   fig-timeseries_hindcast_ENS_acc-v1.py

Run the figure script,

python fig-timeseries_hindcast_ENS_acc-v1.py

The template automatically

reads all ensemble members,
stores the forecast data in

Xens(ensemble, initial, lead)

plots the CTL trajectory,
plots the true state,
plots the three individual ensemble members, and
calculates the ACC using only the deterministic forecast (C00).

The completed figure displays

CTL as a gray line,
True state as black circles,
C00 as orange crosses,
P01 as green crosses, and
M01 as red crosses.

Examine the figure and consider the following questions.

How different are the three ensemble members at the initial time?
Does the distance among the members increase with forecast lead time?
Is one member consistently closer to the true state?
Does the deterministic forecast (C00) always provide the best prediction?

Create the Ensemble-Mean Figure

Next, summarize the ensemble forecasts using the ensemble mean and ensemble spread.

Copy the previous figure,

cp fig-timeseries_hindcast_ENS_acc-v1.py \
   fig-timeseries_hindcast_ENSmean_acc-v1.py

Open

fig-timeseries_hindcast_ENSmean_acc-v1.py

and revise the figure title,

NENS = len(ENSEMBLE_MEMBERS)
figure_title = (
    f"Lorenz Model: Ensemble Mean and Spread (×{NENS})"
)

Calculate the Ensemble Mean

The ensemble forecasts have already been stored in

Xens(ensemble, initial, lead)

Calculate the ensemble mean across the ensemble dimension,

# ============================================================
# Calculate ensemble mean
# ============================================================
Xmean = np.mean(
    Xens,
    axis=0,
)

The resulting array has dimensions

Xmean(initial, lead)

Each element represents the average forecast over the three ensemble members.

Calculate the Ensemble Spread

Next, calculate the ensemble spread,

# ============================================================
# Calculate ensemble spread
# ============================================================
Xstd = np.std(
    Xens,
    axis=0,
)

The resulting array has dimensions

Xstd(initial, lead)

The ensemble spread measures the uncertainty among the ensemble members.

Plot the Ensemble Mean and Spread

Comment out the section that plots the individual ensemble members.

Replace it with

# ============================================================
# Plot ensemble mean and spread
# ============================================================
for iinitial in range(NINITIAL):

    initial_time = float(
        int(HINDCAST_INITIALS[iinitial])
    )

    for ilead in range(NLEAD):

        verification_time = (
            initial_time
            + LEAD_TIMES[ilead]
        )

        if (
            iinitial == 0
            and ilead == 0
        ):
            ensemble_label = (
                "Ensemble mean"
            )
        else:
            ensemble_label = None

        ax[ilead].errorbar(
            verification_time,
            Xmean[
                iinitial,
                ilead,
            ],
            yerr=Xstd[
                iinitial,
                ilead,
            ],
            marker="x",
            markersize=7,
            markeredgewidth=2.0,
            color="C1",
            ecolor="C1",
            elinewidth=1.2,
            capsize=3,
            linestyle="None",
            label=ensemble_label,
        )

The marker represents the ensemble mean, while the error bar represents one standard deviation of the ensemble forecasts.

Calculate ACC Using the Ensemble Mean

The previous figure calculated ACC using only the deterministic forecast (C00).

Replace

ACC[ilead] = np.corrcoef(
    Xtrue[:, ilead],
    Xens[
        iens_c00,
        :,
        ilead,
    ],
)[0, 1]

with

ACC[ilead] = np.corrcoef(
    Xtrue[:, ilead],
    Xmean[:, ilead],
)[0, 1]

Also revise the printed title,

print("")
print("========================================")
print("ACC of Ensemble Mean")
print("========================================")

The remainder of the ACC calculation can be retained without modification.

Run the completed figure,

python fig-timeseries_hindcast_ENSmean_acc-v1.py

Verify the Results

Compare

fig-timeseries_hindcast_ENS_acc-v1.png

and

fig-timeseries_hindcast_ENSmean_acc-v1.png

The first figure displays all individual ensemble members, while the second summarizes the forecasts using the ensemble mean and ensemble spread.

Consider the following questions.

Is the ensemble mean closer to the true state than the deterministic forecast?
Does the ensemble spread generally increase with forecast lead time?
Does the true state remain within the ensemble spread?
Is the ACC of the ensemble mean higher than the ACC of C00?
What additional information does the ensemble spread provide beyond the deterministic forecast?

The ensemble mean provides the best estimate of the forecast state, while the ensemble spread quantifies the uncertainty associated with that forecast. In the following exercises, additional perturbation pairs will be added to investigate how increasing ensemble size influences forecast uncertainty and predictive skill.

Step 5. Increase the Ensemble Size and Evaluate Predictive Skill

In this step, we increase the ensemble size by adding four new ensemble members,

P02
M02
P03
M03

The members P01 and M01 were already generated in Step 3 and therefore do not need to be rerun.

The new perturbation pairs use larger perturbation amplitudes,

P02 and P03    AMP = +1.0
M02 and M03    AMP = -1.0

Together with the existing members,

C00
P01
M01
P02
M02
P03
M03

the complete ensemble contains seven members.

Extend the Ensemble-Hindcast Script

Return to the HINDCAST directory,

cd ../HINDCAST

Copy the ensemble-hindcast script developed in Step 3,

cp run_ENS-v1.sh run_ENS-v2.sh

Open

run_ENS-v2.sh

and revise the ensemble-member settings so that only the four new members are generated:

# ============================================================
# Ensemble-member settings
# ============================================================
MEMBERS=( "P02" "M02" "P03" "M03")

SEEDS=( 2 2 3 3)

AMPS=( 1.0 -1.0 1.0 -1.0)

The member settings are summarized below.

Member	Random seed	Perturbation amplitude
P02	2	+1.0
M02	2	−1.0
P03	3	+1.0
M03	3	−1.0

Because the positive and negative members within each pair use the same random seed and opposite perturbation amplitudes, their initial conditions are symmetric about the observed state.

The remainder of run_ENS-v2.sh can remain unchanged.

Generate the Additional Ensemble Members

Make the script executable,

chmod +x run_ENS-v2.sh

Run the new ensemble members,

./run_ENS-v2.sh

After completion, each initialization directory should contain

HINDCAST/0055/
├── C00/
├── P01/
├── M01/
├── P02/
├── M02/
├── P03/
└── M03/

The same member structure should be present under

Confirm that each newly generated member contains

rest/restart.TTTT.nc
out/

where TTTT is the four-digit initialization time.

Create the Ensemble-Skill Script

Return to the Figure directory,

cd ../Figure

Copy the hindcast-skill template,

cp ../../../Model_src-v3.0/figures/fig-hindcast_skill-template.py \
   fig-hindcast_ENS_skill-v1.py

Open

fig-hindcast_ENS_skill-v1.py

and revise the hindcast-directory template,

FDIR_HINDCAST = ( "../HINDCAST/dddd/member/out/")

Define the Ensemble Members

Add the ensemble-member list,

# ============================================================
# Ensemble members
# ============================================================
ENSEMBLE_MEMBERS = [ "C00", ]

Determine the ensemble size automatically,

NENS = len(ENSEMBLE_MEMBERS)

Revise the figure title using the automatically determined ensemble size,

figure_title = ( f"Lorenz Model: Ensemble Hindcast Skill (x{NENS})")

Also revise the output filename so that the ensemble size is included automatically,

outfile_png = ( f"fig-hindcast_ENS_skill-x{NENS}.png")

For example,

NENS = 1    fig-hindcast_ENS_skill-x1.png
NENS = 3    fig-hindcast_ENS_skill-x3.png
NENS = 5    fig-hindcast_ENS_skill-x5.png
NENS = 7    fig-hindcast_ENS_skill-x7.png

This allows separate figures to be retained for each ensemble size.

Create the Ensemble Forecast Arrays

The original hindcast-skill template stores the forecast as

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

where each experiment corresponds to one hindcast initialization.

To accommodate multiple ensemble members, extend the forecast arrays by adding an ensemble dimension,

Xens(ensemble, experiment, lead)
Yens(ensemble, experiment, lead)
Zens(ensemble, experiment, lead)

The truth arrays remain unchanged,

Xtrue(experiment, lead)
Ytrue(experiment, lead)
Ztrue(experiment, lead)

Define the ensemble dimension,

NENS = len(ENSEMBLE_MEMBERS)

while continuing to use the original

NEXP = len(HINDCAST_EXPERIMENTS)

Create the forecast arrays,

# ============================================================
# Forecast arrays
# ============================================================

Xens = np.full( (NENS, NEXP, NLEAD), np.nan,)
Yens = np.full( (NENS, NEXP, NLEAD), np.nan,)
Zens = np.full( (NENS, NEXP, NLEAD), np.nan,)

Xtrue = np.full( (NEXP, NLEAD), np.nan,)
Ytrue = np.full( (NEXP, NLEAD), np.nan,)
Ztrue = np.full( (NEXP, NLEAD), np.nan,)

The resulting dimensions are

Xens(ensemble, experiment, lead)
Xtrue(experiment, lead)

and similarly for Y and Z.

Read the Ensemble Forecasts

Replace the original hindcast-reading section with nested loops over hindcast experiments and ensemble members.

For each hindcast experiment,

for iexp in range(NEXP):

loop over all ensemble members,

for iens in range(NENS):

and store the forecast into

Xens[iens, iexp, :]

instead of

Xf[iexp, :]

The true states remain identical for all ensemble members and therefore only need to be read once for each hindcast experiment.

Replace

for iexp in range(NEXP):

    experiment = HINDCAST_EXPERIMENTS[iexp]

    ...

    Xf[iexp, :] = Xprd.isel( time=slice(0, NLEAD)).values
    Yf[iexp, :] = Yprd.isel( time=slice(0, NLEAD)).values
    Zf[iexp, :] = Zprd.isel( time=slice(0, NLEAD)).values

with

for iexp in range(NEXP):

    experiment = HINDCAST_EXPERIMENTS[iexp]

    # Convert "0055" to 55.0
    initial_time = float(int(experiment))

    print(
        f"Experiment {experiment}: "
        f"T_INITIAL = {initial_time:.0f}"
    )

    for iens in range(NENS):

        member = ENSEMBLE_MEMBERS[iens]

        FDIR_PRD = FDIR_HINDCAST.replace(
            "dddd",
            experiment,
        ).replace(
            "member",
            member,
        )

        print(
            f"  Read member: {member}"
        )

        fin_prd = xr.open_mfdataset(
            FDIR_PRD + "out.*.nc",
            combine="by_coords",
        )

        time_prd = fin_prd.time
        Xprd = fin_prd.X
        Yprd = fin_prd.Y
        Zprd = fin_prd.Z

        # ----------------------------------------------------
        # Confirm that the hindcast contains enough time steps
        # ----------------------------------------------------

        if time_prd.size < NLEAD:
            raise ValueError(
                f"Experiment {experiment}/{member} "
                f"contains only {time_prd.size} time steps, "
                f"but {NLEAD} time steps are required."
            )

        # ----------------------------------------------------
        # Store ensemble forecast states
        #
        # Xens(ensemble, experiment, lead)
        # ----------------------------------------------------

        Xens[iens, iexp, :] = Xprd.isel( time=slice(0, NLEAD)).values
        Yens[iens, iexp, :] = Yprd.isel( time=slice(0, NLEAD)).values
        Zens[iens, iexp, :] = Zprd.isel( time=slice(0, NLEAD)).values
        fin_prd.close()

After reading all ensemble forecasts, calculate the ensemble mean,

# ============================================================
# Calculate ensemble mean
# ============================================================

Xf = np.mean( Xens, axis=0,)
Yf = np.mean( Yens, axis=0,)
Zf = np.mean( Zens, axis=0,)

The resulting forecast arrays have dimensions

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

which are identical to those used in the original hindcast-skill template.

After calculating the ensemble mean, verify the array dimensions,

print("========================================")
print("Ensemble forecast arrays")
print("========================================")

print("Xens shape :", Xens.shape)
print("Xf shape   :", Xf.shape)
print("Xtrue shape:", Xtrue.shape)

print("")
print(
    "Xens dimensions : "
    "(ensemble, experiment, lead)"
)

print(
    "Xf dimensions   : "
    "(experiment, lead)"
)

print(
    "Xtrue dimensions: "
    "(experiment, lead)"
)

Because Xf, Yf, and Zf have the same dimensions as in the original hindcast-skill template, the remainder of the ACC, RMSE, and RMSE Skill Score calculations can be used without further modification.

Increase the Ensemble Size from One to Three Members

The previous experiment used only the deterministic forecast,

ENSEMBLE_MEMBERS = [ "C00", ]

which produced

NENS = 1

and the figure

fig-hindcast_ENS_skill-x1.png

Next, increase the ensemble size to three members by including the first symmetric perturbation pair,

ENSEMBLE_MEMBERS = [ "C00", "P01", "M01", ]

Because the ensemble size is determined automatically,

NENS = len(ENSEMBLE_MEMBERS)

the program now uses

NENS = 3

and automatically saves the figure as

fig-hindcast_ENS_skill-x3.png

The ensemble mean is still calculated using

Xf = np.mean(Xens, axis=0)
Yf = np.mean(Yens, axis=0)
Zf = np.mean(Zens, axis=0)

However, because Xens, Yens, and Zens now contain three ensemble members, the resulting ensemble mean represents the average of

C00
P01
M01

rather than only the deterministic forecast.

Run the script,

python fig-hindcast_ENS_skill-v1.py

Compare

fig-hindcast_ENS_skill-x1.png

and

fig-hindcast_ENS_skill-x3.png

Consider the following questions.

Does the three-member ensemble improve ACC?
Is RMSE reduced compared with the deterministic forecast?
Does RMSE Skill Score increase?
Which variable (X, Y, or Z) benefits most from ensemble averaging?
Is the improvement similar at all forecast lead times?

The three-member ensemble reduces the influence of individual perturbations by averaging forecasts initialized from slightly different initial conditions. As a result, the ensemble mean often provides a more stable prediction than the deterministic forecast alone.

Key Points

Ensemble forecasting represents forecast uncertainty by integrating multiple forecasts initialized from slightly different initial conditions.
Ensemble forecasts are organized as

Xens(ensemble, experiment, lead)
Yens(ensemble, experiment, lead)
Zens(ensemble, experiment, lead)

The ensemble mean is calculated by averaging over the ensemble dimension,

\[ X_f = \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} X_i, \]

and similarly for \(Y\) and \(Z\).

After averaging, the ensemble-mean forecasts have the same dimensions as the deterministic forecasts developed in Chapter 8,

Xf(experiment, lead)
Yf(experiment, lead)
Zf(experiment, lead)

allowing the existing verification framework to be reused without modification.

The ensemble size is determined automatically from

NENS = len(ENSEMBLE_MEMBERS)

making it straightforward to compare forecasts using different numbers of ensemble members.

Ensemble-mean forecasts are evaluated using the same verification metrics introduced previously:
Anomaly Correlation Coefficient (ACC),
Root Mean Square Error (RMSE), and
RMSE Skill Score (RMSESS).
Increasing the ensemble size generally produces a more stable ensemble mean and often improves predictive skill. However, the improvement typically becomes smaller as additional ensemble members are included.

Homework

In this homework, investigate how increasing the ensemble size influences forecast skill.

Step 1. Generate the Seven-Member Ensemble

Modify

ENSEMBLE_MEMBERS = [ "C00", "P01", "M01", "P02", "M02", "P03", "M03", ]

The ensemble size is determined automatically,

NENS = len(ENSEMBLE_MEMBERS)

which produces

NENS = 7

The ensemble mean is then calculated as

Xf = np.mean(Xens, axis=0)
Yf = np.mean(Yens, axis=0)
Zf = np.mean(Zens, axis=0)

Run

python fig-hindcast_ENS_skill-v1.py

This generates

fig-hindcast_ENS_skill-x7.png

which shows the predictive skill of the seven-member ensemble mean.

Step 2. Compare the Predictive Skill

Compare the seven-member ensemble results with the deterministic forecast (x1) developed in this chapter.

Consider the following questions.

Does the seven-member ensemble improve ACC?
Does RMSE decrease?
Does RMSE Skill Score increase?
Which variable (X, Y, or Z) benefits most from increasing the ensemble size?
Is the improvement similar at all forecast lead times?

Submission

Submit

fig-hindcast_ENS_skill-x7.png
A short paragraph (approximately 100–200 words) comparing the predictive skill of the seven-member ensemble with the one-member (C00) forecast.

Your discussion should explain how increasing the ensemble size influences the forecast skill and identify which verification metrics show the largest improvement.