Chapter 9.1: Ensemble Forecasting I: Perturbed Initial Conditions

Learning Objectives

After completing this chapter, students will be able to:

Explain why observational uncertainty leads to uncertainty in the initial conditions of a forecast model.
Describe the rationale for ensemble forecasting using perturbed initial conditions.
Generate symmetric random perturbations around an observation to construct an ensemble of initial conditions.
Perform ensemble forecast experiments using the Lorenz climate model.
Explain how the perturbation amplitude influences the range of possible forecast trajectories and forecast uncertainty.
Recognize common operational approaches for generating ensemble initial conditions, including the Breeding Method, Singular Vectors (SVs), Ensemble Kalman Filter (EnKF), and Lagged Average Forecasting (LAF).

Concept

Why do we need perturbed initial conditions?

In the previous chapter, we initialized the forecast model directly from a single observation and performed a deterministic hindcast experiment. However, observations inevitably contain measurement errors and sampling uncertainty, meaning that the observed state is only an estimate of the unknown true state of the climate system. Consequently, a deterministic forecast initialized from a single observation represents only one possible future evolution.

Ensemble forecasting addresses this problem by generating multiple initial conditions that are consistent with the observational uncertainty. Instead of assuming that the observation is exact, the model is initialized from a collection of slightly different initial conditions, producing multiple forecast realizations that collectively represent a range of plausible future states.

Operational forecasting centers employ sophisticated perturbation techniques, such as the Breeding Method (Bred Vectors), Singular Vectors (SVs), and the Ensemble Kalman Filter (EnKF), to generate dynamically consistent perturbations that efficiently represent uncertainty in the atmospheric or climate state. In this course, however, we adopt a much simpler approach by generating symmetric random perturbations around the observation. Although this approach is less sophisticated than operational techniques, it captures the fundamental concept of ensemble forecasting while remaining computationally simple and easy to implement.

To preserve the observation as the best estimate of the initial state, perturbations are generated in positive–negative pairs (e.g., +δ and −δ). This symmetric perturbation strategy ensures that the ensemble mean is identical to the observation at the initial time while representing the uncertainty surrounding that estimate. Throughout this chapter, the unperturbed forecast is referred to as the control (CTL) experiment, while the positively and negatively perturbed members are denoted as P01, M01, P02, M02, and so forth.

An important consideration is the perturbation amplitude. If the perturbations are too small, the ensemble members remain tightly clustered and may fail to encompass the true initial state, resulting in an overconfident forecast. Conversely, excessively large perturbations increase the likelihood of including the true state but also introduce unrealistic initial conditions and unnecessarily large forecast uncertainty. Therefore, the perturbation amplitude should be selected to realistically represent the expected observational uncertainty.

Alternative approach: Lagged Average Forecasting (LAF)

Perturbing the initial conditions is not the only way to construct an ensemble forecast. Another widely used approach is Lagged Average Forecasting (LAF), also known as a lagged ensemble. Instead of generating multiple perturbed initial states at a single initialization time, the lagged ensemble combines forecasts initialized at different times (e.g., forecasts started 6, 12, or 24 hours earlier) into a single ensemble.

Because these forecasts are already available from routine operational model runs, the LAF approach can substantially increase the ensemble size with little additional computational cost. However, lagged ensemble members originate from different analyses and forecast lead times, so they do not represent initial-condition uncertainty as directly as perturbation-based ensembles. In practice, many operational numerical weather prediction and climate prediction systems combine multiple ensemble generation strategies to achieve an effective balance between forecast quality and computational efficiency.

Model/Equation

Observation Model

Following Chapter 8, the synthetic observation is assumed to contain random measurement errors,

\[ \mathbf{y} = \mathbf{x} + \boldsymbol{\varepsilon}, \]

where

\(\mathbf{x}\) is the unknown true state,
\(\mathbf{y}\) is the observation, and
\(\boldsymbol{\varepsilon}\) represents the observation error.

Because the true state is unknown, the observation serves as the best available estimate of the initial condition for the forecast experiment.

Ensemble Initialization

The control simulation (CTL) generated in the previous chapters continues to represent the true climate evolution and serves as the reference solution for forecast verification. The ensemble forecasts are initialized from the observation rather than the true state.

The central ensemble member (C00) is initialized directly from the observation,

\[ \mathbf{x}_{\mathrm{C00}} = \mathbf{y}. \]

Additional ensemble members are generated by adding and subtracting symmetric random perturbations,

\[ \mathbf{x}_{\mathrm{P}i} = \mathbf{y} + \boldsymbol{\delta}_i, \]

\[ \mathbf{x}_{\mathrm{M}i} = \mathbf{y} - \boldsymbol{\delta}_i, \]

where

\(\boldsymbol{\delta}_i\) is the \(i\)-th random perturbation,
P\(i\) denotes the positively perturbed member, and
M\(i\) denotes the negatively perturbed member.

The perturbations are generated in positive–negative pairs so that

\[ \sum_{i=1}^{N} \boldsymbol{\delta}_i - \sum_{i=1}^{N} \boldsymbol{\delta}_i = \mathbf{0}, \]

which guarantees that the ensemble mean at the initial time is equal to the observation,

\[ \overline{\mathbf{x}}(t_0) = \mathbf{y}. \]

Thus, the observation remains the best estimate of the initial condition while the ensemble members represent the uncertainty surrounding that estimate.

Random Perturbation Model

In this chapter, each perturbation is independently sampled from a Gaussian distribution,

\[ \boldsymbol{\delta} \sim \mathcal{N} \left( \mathbf{0}, \sigma_p^2 \mathbf{I} \right), \]

where

\(\sigma_p\) is the prescribed perturbation amplitude, and
\(\mathbf{I}\) is the identity matrix.

The perturbation amplitude controls the initial spread of the ensemble. Smaller perturbations produce a narrow ensemble with high confidence but may fail to encompass the true initial state. Larger perturbations increase the likelihood of including the true state but also produce a wider ensemble and greater forecast uncertainty.

Ensemble Member Naming Convention

Throughout this chapter, the ensemble members are named according to the following convention.

Member	Description	Initial Condition
CTL	Perfect-model (truth) simulation	True state
C00	Central ensemble member	\(\mathbf{y}\)
P01	Positive perturbation #1	\(\mathbf{y}+\boldsymbol{\delta}_1\)
M01	Negative perturbation #1	\(\mathbf{y}-\boldsymbol{\delta}_1\)
P02	Positive perturbation #2	\(\mathbf{y}+\boldsymbol{\delta}_2\)
M02	Negative perturbation #2	\(\mathbf{y}-\boldsymbol{\delta}_2\)
\(\vdots\)	\(\vdots\)	\(\vdots\)

This symmetric perturbation strategy preserves the observation as the ensemble mean while generating multiple plausible initial conditions for the forecast experiment.

Sample Code

The sample programs in this section reuse the Model_src-v3.0 framework developed in the previous exercises. The control simulation provides the true climate state, while the ensemble forecasts are initialized from synthetic observations with prescribed observation error.

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_obs.py
│   └── fig-timeseries_CTL_PRD-template.py
├── utils/
│   ├── get_observation.py
│   └── create_initial_obs.py
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

The individual source files used in this section are available 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.
`get_observation.py`	Generates synthetic observations from the control simulation by adding random observation error.
`create_initial_obs.py`	Creates a forecast initial condition directly from the synthetic observations at a prescribed initialization time.
`fig-timeseries_obs.py`	Plots the control simulation together with the synthetic observations.
`fig-timeseries_CTL_PRD-template.py`	Provides a template for comparing the control simulation with a forecast initialized from the observations.

The control simulation is integrated over

TSTR = 0
TEND = 100

and synthetic observations are generated over

OBS_TSTR = 55
OBS_TEND = 95

using the prescribed observation error

OERR = 1.0

The ensemble forecast case study is initialized at

T_INITIAL = 60

and integrated over

TSTR = 60
TEND = 65

The ensemble experiments are stored in separate directories according to their member names:

EXP05/
├── C00/
├── P01/
├── M01/
├── P02/
├── M02/
├── P03/
├── M03/
└── run_member.sh

The central member (C00) is initialized directly from the synthetic observation. The remaining ensemble members are generated by adding symmetric random perturbations to the observed initial condition:

P01 and M01    seed = 1, amplitude = ±0.5
P02 and M02    seed = 2, amplitude = ±0.5
P03 and M03    seed = 3, amplitude = ±0.5

The positively and negatively perturbed members within each pair use the same random seed and opposite perturbation amplitudes. Therefore, the perturbations cancel when averaged across each pair, and the ensemble mean at the initial time remains equal to the observed state.

During the exercises, the original source programs are progressively extended to create

create_initial_obs-perturbed.py
run_member.sh
fig-timeseries_CTL_C00.py
fig-timeseries_CTL_ENS.py
fig-timeseries_CTL_ENSmean.py

These programs are developed as part of the exercises and are therefore not included in the original Model_src-v3.0 source-code archive.

The individual forecast members are reorganized into two-dimensional arrays such as

Xens(ensemble, time)
Yens(ensemble, time)
Zens(ensemble, time)

The first dimension represents the ensemble member, while the second dimension represents forecast time. This data structure makes it straightforward to calculate ensemble statistics across all forecast members.

The following two ensemble statistics are evaluated:

Ensemble Mean represents the average forecast trajectory across all ensemble members.
Ensemble Spread measures the uncertainty of the ensemble using the standard deviation of the ensemble members.

The ensemble spread is interpreted as

Small spread    high agreement among ensemble members
Large spread    large forecast uncertainty

The final ensemble figure displays the control simulation together with the ensemble mean and the ensemble spread (± one standard deviation), illustrating both the expected forecast evolution and its associated uncertainty.

Exercise

Step 1. Generate the Control Simulation and Synthetic Observations

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

Create the working directories,

mkdir Chapter9_1
cd Chapter9_1

mkdir CTL
cd CTL
cp ../../Model_src-v3.0/*.* .

The CTL 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 generates the reference simulation over the period \(t=0\)–100.

Generate Synthetic Observations

Synthetic observations are generated by sampling the CTL simulation and adding random observation errors. In this chapter, use an observation error of

OERR = 1.0

to clearly illustrate uncertainty in the observed initial conditions.

Copy the observation-generation script into the CTL directory,

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

Open get_observation.py and configure the analysis period and observation period as

TSTR = 50
TEND = 100

OBS_TSTR = 55.0
OBS_TEND = 95.0

OERR = 1.0

Then execute

python get_observation.py

This creates

observation.nc

which contains synthetic observations that will be used to initialize the ensemble forecast experiments.

Visualize the CTL Simulation and Observations

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

cd ../

mkdir Figure
cd Figure

Copy the observation time-series script,

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

Open

fig-timeseries_obs.py

and configure the model output directory, observation file, and plotting period as

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

TSTR = 50
TEND = 95

Run the figure script,

python fig-timeseries_obs.py

The figure should display

the CTL simulation, representing the true evolution of the Lorenz system, and
the synthetic observations, which differ from the true state because of the prescribed observation error.

Verify the Results

After completing this step, the directory structure should resemble

Chapter9_1/
├── 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 uncertain observations that will be used to construct the ensemble initial conditions.

Step 2. Generate the Central Ensemble Member (C00)

Before generating perturbed ensemble members, we first initialize the model directly from an observation. This experiment, referred to as C00, represents the central ensemble member and provides the unperturbed forecast around which the other ensemble members will be constructed.

Return to the Chapter9_1 directory,

cd ../

Create the experiment directory,

mkdir EXP05
cd EXP05

Copy the initial-condition generator into EXP05,

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

Next, create the directory for the central ensemble member,

mkdir C00
cd C00

Copy the required model source files,

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

Copy the initial-condition generator from the EXP05 directory,

cp ../create_initial_obs.py .

Generate the Initial Condition

Generate the restart file using the observation at \(t=60\),

python create_initial_obs.py  --t_initial=60

This creates

rest/restart.0060.nc

The restart file is initialized directly from the synthetic observation without an additional perturbation.

Run the C00 Forecast

Open run_model_climate-v3.sh and configure the forecast as

REST=True
TSTR=60
TEND=65

Then execute

bash run_model_climate-v3.sh

This deterministic forecast serves as the central ensemble member (C00).

Visualize and Verify the C00 Forecast

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

cd ../../

mkdir Figure
cd Figure

Copy the time-series comparison template,

cp ../../../Model_src-v3.0/figures/fig-timeseries_CTL_PRD-template.py \
   fig-timeseries_CTL_C00.py

Open

fig-timeseries_CTL_C00.py

and configure the model output directories as

# Model output directories
# Model output directories
FDIR_CTL = "../CTL/out/"
FDIR_PRD = "../EXP05/C00/out/"

Use a longer period for the CTL simulation and the forecast period for C00,

# Control simulation period
TSTR_CTL = 55
TEND_CTL = 70

# Forecast period
TSTR_PRD = 60
TEND_PRD = 65

Revise the figure title and legend so that the forecast is labeled C00 rather than PRD.

Run the figure script,

python fig-timeseries_CTL_C00.py

The figure compares

the CTL simulation, representing the true evolution of the Lorenz system, and
the C00 forecast, initialized directly from the observation at \(t=60\).

Examine the initial difference between CTL and C00 and how this difference evolves over the forecast period. Because the synthetic observation contains random error, C00 does not necessarily begin from the true state.

Verify the Results

After completing this step, the directory structure should resemble

Chapter9_1/
├── CTL/
│   ├── out/
│   ├── rest/
│   ├── observation.nc
│   ├── get_observation.py
│   ├── Lmodel_climate-v3.py
│   └── run_model_climate-v3.sh
├── Figure/
│   ├── fig-timeseries_obs.py
    └── fig-timeseries_CTL_C00.py
└── EXP05/
    ├── C00/
    │   ├── out/
    │   ├── rest/
    │   │   └── restart.0070.nc
    │   ├── create_initial_obs.py
    │   ├── Lmodel_climate-v3.py
    │   └── run_model_climate-v3.sh
    └── create_initial_obs.py

In the next step, we will modify the initial-condition generator to add a random perturbation and create the first perturbed ensemble member, P01.

Step 3. Generate the First Perturbed Member (P01)

In this step, we modify the initial-condition generator to add a random perturbation to the observed initial state. The modified script is then used to generate the first positively perturbed ensemble member, P01.

Return to the EXP05 directory,

cd ../EXP05

Copy the original initial-condition generator to a new file,

cp create_initial_obs.py create_initial_obs-perturbed.py

The original create_initial_obs.py remains unchanged and is used for the central ensemble member (C00). The new script, create_initial_obs-perturbed.py, will be used to generate all perturbed ensemble members.

Modify the Initial-Condition Generator

Open

create_initial_obs-perturbed.py

and add two new command-line arguments,

--seed: random number seed
--amp: perturbation amplitude

For example,

parser.add_argument(
    "--seed",
    type=int,
    default=1,
    help="Random seed",
)

parser.add_argument(
    "--amp",
    type=float,
    default=0.0,
    help="Perturbation amplitude",
)

After parsing the command-line arguments, define

SEED = args.seed
AMP = args.amp

Initialize the random number generator,

rng = np.random.default_rng(SEED)

and generate three independent Gaussian perturbations,

dX = rng.normal(0.0, 1.0)
dY = rng.normal(0.0, 1.0)
dZ = rng.normal(0.0, 1.0)

Replace the original observed initial conditions with

Xobs_value = obs.Xobs.item()
Yobs_value = obs.Yobs.item()
Zobs_value = obs.Zobs.item()

Xint_value = Xobs_value + AMP * dX
Yint_value = Yobs_value + AMP * dY
Zint_value = Zobs_value + AMP * dZ

The random seed determines the direction of the perturbation, while AMP controls its magnitude and sign.

You may also add the random seed and perturbation amplitude to the restart-file attributes and printed summary.

Generate the P01 Forecast

Create the directory,

mkdir P01
cd P01

Copy the required files,

cp ../C00/Lmodel_climate-v3.py .
cp ../C00/run_model_climate-v3.sh .
cp ../create_initial_obs-perturbed.py .

Generate the perturbed restart file,

python create_initial_obs-perturbed.py \
    --t_initial=60 \
    --seed=1 \
    --amp=0.5

This creates

rest/restart.0060.nc

for the P01 ensemble member.

Open run_model_climate-v3.sh and configure

REST=True
TSTR=60
TEND=65

Then execute

bash run_model_climate-v3.sh

Visualize the Forecast Ensemble

Return to the Figure directory,

cd ../../Figure

Copy the C00 figure script,

cp fig-timeseries_CTL_C00.py \
   fig-timeseries_CTL_ENS.py

Open

fig-timeseries_CTL_ENS.py

and modify the forecast directories,

FDIR_CTL = "../CTL/out/"

FDIR_ENS = {
    "C00": "../EXP05/C00/out/",
    "P01": "../EXP05/P01/out/",
}

Use

TSTR_CTL = 55
TEND_CTL = 70

TSTR_ENS = 60
TEND_ENS = 65

Read the Ensemble Members

Replace the section that reads the single perturbed forecast (Xprd, Yprd, and Zprd) with the following code:

# ============================================================
# Read ensemble simulations
# ============================================================
ensemble_members = list(FDIR_ENS.keys())

Xens_list = []
Yens_list = []
Zens_list = []

for member in ensemble_members:

    fdir = FDIR_ENS[member]

    fin_ens = xr.open_mfdataset(
        fdir + "out.*.nc",
        combine="by_coords",
    )

    Xtmp = fin_ens.X.sel(
        time=slice(TSTR_ENS, TEND_ENS)
    ).load()

    Ytmp = fin_ens.Y.sel(
        time=slice(TSTR_ENS, TEND_ENS)
    ).load()

    Ztmp = fin_ens.Z.sel(
        time=slice(TSTR_ENS, TEND_ENS)
    ).load()

    Xens_list.append(Xtmp)
    Yens_list.append(Ytmp)
    Zens_list.append(Ztmp)

    fin_ens.close()


# ============================================================
# Combine ensemble members
# ============================================================
Xens = xr.concat(
    Xens_list,
    dim=xr.IndexVariable(
        "ensemble",
        ensemble_members,
    ),
)

Yens = xr.concat(
    Yens_list,
    dim=xr.IndexVariable(
        "ensemble",
        ensemble_members,
    ),
)

Zens = xr.concat(
    Zens_list,
    dim=xr.IndexVariable(
        "ensemble",
        ensemble_members,
    ),
)

time_ens = Xens.time


# ============================================================
# Print ensemble information
# ============================================================
print("Xens dimensions:", Xens.dims)
print("Xens shape     :", Xens.shape)
print("Members        :", Xens.ensemble.values)

The resulting ensemble arrays have the dimensions

Xens(ensemble, time)
Yens(ensemble, time)
Zens(ensemble, time)

For this step, the ensemble coordinate contains

C00
P01

Because each dataset is closed inside the loop, remove or comment out the original closing command for the perturbed simulation,

# fin_prd.close()

The control dataset can still be closed at the end of the script,

fin_ctl.close()

Plot the Ensemble Members

Replace the plotting command for the single perturbed forecast with a loop over the ensemble members. Use black for the CTL simulation and the default Matplotlib color cycle for the forecast members.

For the X component,

ax[0].plot(
    time_ctl,
    Xctl,
    color="black",
    linewidth=2.0,
    label="CTL",
)

for member in ensemble_members:

    ax[0].plot(
        time_ens,
        Xens.sel(ensemble=member),
        linewidth=1.0,
        label=member,
    )

Apply the same modification to the Y and Z panels.

Revise the figure title as

ax[0].set_title(
    "Control Simulation and Ensemble Forecast"
)

The completed figure should display

CTL: true evolution of the Lorenz system,
C00: forecast initialized directly from the observation,
P01: forecast initialized from the positively perturbed observation.

Run

python fig-timeseries_CTL_ENS.py

Verify the Results

Compare the trajectories of CTL, C00, and P01.

Consider the following questions.

How different are the initial conditions of C00 and P01?
Does the difference between C00 and P01 increase during the forecast?
Is the P01 trajectory closer to or farther from the CTL simulation than C00?
Can two forecast members adequately represent the uncertainty in the observed initial condition?

At this stage, the forecast ensemble contains only

C00
P01

and therefore is not yet symmetric about the observation. In the next step, we generate the corresponding negatively perturbed member (M01) using the same random perturbation with the opposite sign.

Step 4. Automate the Generation of the M01 Member

Rather than manually repeating the same commands, we develop a shell script that automatically generates one ensemble member. By modifying only a few parameters, the same script can later be reused to generate additional ensemble members.

Return to the EXP05 directory,

cd ../EXP05

Create a shell script,

run_member.sh

and add the following code.

#!/bin/bash

MEMBER="M01"

T_INITIAL=60
T_END=65

SEED=1
AMP=-0.5

mkdir -p ${MEMBER}
cd ${MEMBER}

cp ../C00/Lmodel_climate-v3.py .
cp ../C00/run_model_climate-v3.sh .
cp ../create_initial_obs-perturbed.py .

python create_initial_obs-perturbed.py \
    --t_initial=${T_INITIAL} \
    --seed=${SEED} \
    --amp=${AMP}

bash run_model_climate-v3.sh

cd ../

Make the script executable,

chmod +x run_member.sh

Run

./run_member.sh

Visualize the Three-Member Ensemble

Return to the Figure directory,

cd ../Figure

Open

fig-timeseries_CTL_ENS.py

and update the ensemble member list,

FDIR_ENS = {
    "C00": "../EXP05/C00/out/",
    "P01": "../EXP05/P01/out/",
    "M01": "../EXP05/M01/out/",
}

Run

python fig-timeseries_CTL_ENS.py

Verify the Results

Compare the trajectories of CTL, C00, P01, and M01.

Consider the following questions.

Are P01 and M01 approximately symmetric about C00 at the initial time?
Does the ensemble spread increase with forecast lead time?
Does the three-member ensemble better represent uncertainty than the two-member ensemble?

After completing this step, the directory structure should resemble

Chapter9_1/
├── CTL/
├── Figure/
│   ├── fig-timeseries_obs.py
│   ├── fig-timeseries_CTL_C00.py
│   └── fig-timeseries_CTL_ENS.py
└── EXP05/
    ├── C00/
    ├── P01/
    ├── M01/
    ├── create_initial_obs.py
    ├── create_initial_obs-perturbed.py
    ├── Lmodel_climate-v3.py
    ├── run_model_climate-v3.sh
    └── run_member.sh

The forecast ensemble now consists of the central member (C00) and one symmetric perturbation pair (P01 and M01). In Chapter 9.2, this workflow will be extended to generate larger ensembles and investigate how ensemble size influences forecast uncertainty and predictive skill.

Step 5. Generate Additional Ensemble Members and Visualize the Ensemble Mean and Spread

In this step, we reuse run_member.sh to generate four additional ensemble members,

P02
M02
P03
M03

Together with C00, P01, and M01, these experiments form a seven-member forecast ensemble,

C00
P01  M01
P02  M02
P03  M03

Each positive–negative pair uses the same random seed and opposite perturbation amplitudes. Therefore, the initial perturbations are symmetric about the observed initial condition.

Generate Additional Ensemble Members

Return to the EXP05 directory,

cd ../EXP05

Open

run_member.sh

and change the ensemble-member settings for each experiment.

Generate P02 using

MEMBER="P02"

T_INITIAL=60
T_END=65

SEED=2
AMP=0.5

Then execute

./run_member.sh

Next, generate M02 using the same random seed and the opposite perturbation amplitude,

MEMBER="M02"

T_INITIAL=60
T_END=65

SEED=2
AMP=-0.5

Run

./run_member.sh

Generate P03 using

MEMBER="P03"

T_INITIAL=60
T_END=65

SEED=3
AMP=0.5

Run

./run_member.sh

Finally, generate M03 using

MEMBER="M03"

T_INITIAL=60
T_END=65

SEED=3
AMP=-0.5

Run

./run_member.sh

The ensemble-member settings are summarized below.

Member	Random seed	Perturbation amplitude
C00	—	0.0
P01	1	+0.5
M01	1	−0.5
P02	2	+0.5
M02	2	−0.5
P03	3	+0.5
M03	3	−0.5

Because each positive–negative pair uses the same random seed, the two members in each pair use the same random perturbation vector with opposite signs.

Verify the Forecast Experiments

Confirm that the following directories were created under EXP05:

EXP05/
├── C00/
├── P01/
├── M01/
├── P02/
├── M02/
├── P03/
├── M03/
├── create_initial_obs.py
├── create_initial_obs-perturbed.py
└── run_member.sh

Each newly generated ensemble-member directory should contain

rest/restart.0060.nc
out/

The complete ensemble now contains seven members,

\[ N_{\mathrm{ens}} = 7. \]

At the initial time, the positive and negative perturbations cancel within each pair. Therefore, the seven-member ensemble mean is equal to the observed initial condition.

Create the Ensemble-Mean Figure

Return to the Figure directory,

cd ../Figure

Copy the ensemble time-series script,

cp fig-timeseries_CTL_ENS.py \
   fig-timeseries_CTL_ENSmean.py

Open

fig-timeseries_CTL_ENSmean.py

and update the ensemble-member directories,

FDIR_CTL = "../CTL/out/"

FDIR_ENS = {
    "C00": "../EXP05/C00/out/",
    "P01": "../EXP05/P01/out/",
    "M01": "../EXP05/M01/out/",
    "P02": "../EXP05/P02/out/",
    "M02": "../EXP05/M02/out/",
    "P03": "../EXP05/P03/out/",
    "M03": "../EXP05/M03/out/",
}

Continue to use

TSTR_CTL = 55
TEND_CTL = 70

TSTR_ENS = 60
TEND_ENS = 65

The existing code reads the individual forecast members and creates the following arrays:

Xens(ensemble, time)
Yens(ensemble, time)
Zens(ensemble, time)

Calculate the Ensemble Mean and Spread

After creating Xens, Yens, and Zens, calculate the ensemble mean,

# ============================================================
# Calculate ensemble mean
# ============================================================
Xmean = Xens.mean(dim="ensemble")
Ymean = Yens.mean(dim="ensemble")
Zmean = Zens.mean(dim="ensemble")

Calculate the ensemble spread using the standard deviation across the ensemble members,

# ============================================================
# Calculate ensemble spread
# ============================================================
Xstd = Xens.std(dim="ensemble")
Ystd = Yens.std(dim="ensemble")
Zstd = Zens.std(dim="ensemble")

The ensemble mean represents the average forecast trajectory, while the ensemble standard deviation measures the spread of the forecast members around the ensemble mean.

For example,

\[ \overline{X}(t) = \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} X_i(t), \]

and

\[ \sigma_X(t) = \sqrt{ \frac{1}{N_{\mathrm{ens}}} \sum_{i=1}^{N_{\mathrm{ens}}} \left[ X_i(t)-\overline{X}(t) \right]^2 }. \]

The same calculations are applied to \(Y\) and \(Z\).

Add the Ensemble Mean and Spread

First, retain the individual ensemble-member trajectories and add the ensemble mean and spread to the existing figure.

For the \(X\) component, add

ax[0].fill_between(
    time_ens,
    Xmean - Xstd,
    Xmean + Xstd,
    color="C1",
    alpha=0.2,
    label="Ensemble spread",
)

ax[0].plot(
    time_ens,
    Xmean,
    color="C1",
    linewidth=2.0,
    label="Ensemble mean",
)

Run the figure script,

python fig-timeseries_CTL_ENSmean.py

Examine how the ensemble mean and shaded spread summarize the individual forecast trajectories.

The shaded region represents

\[ \overline{X}(t) \pm \sigma_X(t), \]

rather than the complete minimum-to-maximum range of the ensemble.

Simplify the Ensemble Figure

Next, remove or comment out the plotting commands for the individual ensemble members. Retain only the CTL simulation, ensemble mean, and ensemble spread.

For example,

# ============================================================
# X component
# ============================================================
# for member in ensemble_members:
#     ax[0].plot(
#         time_ens,
#         Xens.sel(ensemble=member),
#         linewidth=1.0,
#         label=member,
#     )

ax[0].fill_between(
    time_ens,
    Xmean - Xstd,
    Xmean + Xstd,
    color="C1",
    alpha=0.2,
    label="Ensemble spread",
)

ax[0].plot(
    time_ens,
    Xmean,
    color="C1",
    linewidth=2.0,
    label="Ensemble mean",
)

Apply the same procedure to the \(Y\) component,

# ============================================================
# Y component
# ============================================================
# for member in ensemble_members:
#     ax[1].plot(
#         time_ens,
#         Yens.sel(ensemble=member),
#         linewidth=1.0,
#         label=member,
#     )

ax[1].fill_between(
    time_ens,
    Ymean - Ystd,
    Ymean + Ystd,
    color="C1",
    alpha=0.2,
)

ax[1].plot(
    time_ens,
    Ymean,
    color="C1",
    linewidth=2.0,
    label="Ensemble mean",
)

Apply the same procedure to the \(Z\) component,

# ============================================================
# Z component
# ============================================================
# for member in ensemble_members:
#     ax[2].plot(
#         time_ens,
#         Zens.sel(ensemble=member),
#         linewidth=1.0,
#         label=member,
#     )

ax[2].fill_between(
    time_ens,
    Zmean - Zstd,
    Zmean + Zstd,
    color="C1",
    alpha=0.2,
)

ax[2].plot(
    time_ens,
    Zmean,
    color="C1",
    linewidth=2.0,
    label="Ensemble mean",
)

The final figure should display only

CTL
Ensemble spread
Ensemble mean

Run the revised figure script,

python fig-timeseries_CTL_ENSmean.py

Verify the Results

The completed figure should display

the CTL trajectory,
the ensemble mean, and
the ensemble mean ± one standard deviation as the shaded ensemble spread.

Consider the following questions.

Is the ensemble mean equal to the observed initial condition at \(t=60\)?
Does the ensemble spread increase with forecast lead time?
Does the CTL trajectory remain within the shaded ensemble spread?
Is the ensemble mean closer to the CTL trajectory than C00?
Are there periods when the ensemble spread is narrow but the ensemble mean is far from CTL?
What do these results suggest about forecast uncertainty and predictability?

The ensemble spread represents the degree of agreement among the ensemble members. A narrow spread indicates that the members produce similar forecasts, whereas a wide spread indicates greater forecast uncertainty. However, a narrow spread does not necessarily indicate an accurate forecast if the complete ensemble misses the true trajectory.

After completing this step, the directory structure should resemble

Chapter9_1/
├── CTL/
├── Figure/
│   ├── fig-timeseries_obs.py
│   ├── fig-timeseries_CTL_C00.py
│   ├── fig-timeseries_CTL_ENS.py
│   └── fig-timeseries_CTL_ENSmean.py
└── EXP05/
    ├── C00/
    ├── P01/
    ├── M01/
    ├── P02/
    ├── M02/
    ├── P03/
    ├── M03/
    ├── create_initial_obs.py
    ├── create_initial_obs-perturbed.py
    └── run_member.sh

The seven-member ensemble provides a broader representation of the uncertainty surrounding the observed initial condition than the three-member ensemble. In the next chapter, this workflow will be extended to investigate how the number of ensemble members influences forecast stability and predictive skill.

Key Points

Ensemble forecasting represents uncertainty in the observed initial condition by integrating multiple forecasts from slightly different initial conditions.
In this chapter, ensemble members are generated by adding random Gaussian perturbations to the synthetic observations.
The central ensemble member (C00) is initialized directly from the observation, while the remaining members are created by adding positive and negative perturbations.
Each positive–negative perturbation pair uses the same random seed with opposite perturbation amplitudes. As a result, the perturbations are symmetric about the observed initial condition, and the ensemble mean at the initial time remains equal to the observation.
The perturbation amplitude controls the initial ensemble spread. A larger perturbation amplitude produces a wider range of initial conditions and generally leads to larger forecast uncertainty.
Organizing the forecast trajectories into arrays such as

Xens(ensemble, time)
Yens(ensemble, time)
Zens(ensemble, time)

makes it straightforward to calculate ensemble statistics.

The ensemble mean is calculated by averaging all ensemble members at each forecast time.
The ensemble spread is measured by the standard deviation of the ensemble members and quantifies the uncertainty of the ensemble forecast.
The ensemble mean provides a single best estimate of the forecast evolution, while the ensemble spread indicates the degree of agreement among the forecast members.
A small ensemble spread indicates that the forecast members are similar, whereas a large ensemble spread indicates greater forecast uncertainty. However, a small spread does not necessarily imply an accurate forecast if the entire ensemble misses the true trajectory.
In this chapter, only random perturbations are used to generate the ensemble members. More sophisticated methods, such as Lagged Average Forecasting (LAF), Breeding Vectors (BV), Singular Vectors (SV), and the Ensemble Kalman Filter (EnKF), will be introduced in later discussions.

Homework

Repeat the ensemble forecast experiment using a larger perturbation amplitude.

Step 1. Create a New Experiment

Create a new experiment directory,

mkdir EXP20
cd EXP20

Copy the central ensemble member (C00) and the required source files,

cp -r ../EXP05/C00 .
cp ../EXP05/create_initial_obs.py .
cp ../EXP05/create_initial_obs-perturbed.py .
cp ../EXP05/run_member.sh .

The directory should resemble

EXP20/
├── C00/
├── create_initial_obs.py
├── create_initial_obs-perturbed.py
└── run_member.sh

Step 2. Increase the Perturbation Amplitude

Modify run_member.sh so that every perturbed ensemble member uses

AMP=2.0

for the positive members and

AMP=-2.0

for the negative members.

Retain the same initialization time,

T_INITIAL = 60

and use the same random seeds developed in this chapter.

Step 3. Generate Additional Ensemble Members

Generate the six perturbed ensemble members,

P01
M01
P02
M02
P03
M03

using run_member.sh.

The completed directory should resemble

EXP20/
├── C00/
├── P01/
├── M01/
├── P02/
├── M02/
├── P03/
├── M03/
├── create_initial_obs.py
├── create_initial_obs-perturbed.py
└── run_member.sh

Step 4. Calculate the Ensemble Mean and Spread

Copy the ensemble-mean figure script,

cp ../Figure/fig-timeseries_CTL_ENSmean.py \
   ../Figure/fig-timeseries_CTL_ENSmean_EXP20.py

Modify the ensemble directories so that the script reads the ensemble members under EXP20, then generate the ensemble mean and spread.

Questions

Compare the results of EXP05 and EXP20.

Discuss the following questions.

How does increasing the perturbation amplitude affect the initial ensemble spread?
Does the ensemble spread grow more rapidly during the forecast?
Does the ensemble mean remain similar between the two experiments?
Does a larger initial perturbation necessarily produce a better forecast?

Submission

Submit

Figure: fig-timeseries_CTL_ENSmean_EXP20.png.
Discussion: A brief comparison (approximately 150–250 words) describing the differences in ensemble spread and ensemble mean between EXP05 and EXP20, and explaining how increasing the perturbation amplitude influences forecast uncertainty.