8.1 Initial Error Growth and Predictability

Learning Objectives

After completing this lecture, students should be able to:

explain why weather forecasts have a finite predictability limit,
distinguish between error growth in linear and nonlinear systems,
reproduce the classic Lorenz initial-condition sensitivity experiment,
quantify forecast error using two simulations with slightly different initial conditions, and
interpret how chaotic dynamics amplify tiny initial-condition errors.

Concept

The previous chapter demonstrated that different initial conditions produce different climate realizations in an ensemble simulation. Although every ensemble member used the same governing equations and the same model configuration, their trajectories gradually diverged because each simulation began from a slightly different initial state.

This sensitivity to initial conditions also determines the predictability of weather forecasts. Since atmospheric observations are never perfectly accurate, every numerical weather prediction begins with small uncertainties in the initial state of the atmosphere. The important question is how these small errors evolve as the forecast progresses.

To understand this process, we first examine a simple linear model in which a small initial perturbation produces a predictable response. We then repeat the same experiment using the Lorenz climate model introduced in Section 7.1.

The Lorenz model became famous through the work of Edward Lorenz in 1961. While repeating one of his numerical experiments, Lorenz restarted the model using values rounded to fewer decimal places than those stored internally by the computer. Although the difference between the two initial conditions was extremely small, the simulations eventually evolved into completely different trajectories.

This experiment demonstrated that deterministic systems can nevertheless possess limited predictability. Even with a perfect numerical model, weather forecasts cannot remain accurate indefinitely because the atmosphere's initial state can never be observed perfectly.

In this section, we reproduce Lorenz's experiment using restart files. After completing a control simulation, we create a second restart file by modifying only the least significant digits of the model state. Comparing the original and perturbed simulations illustrates how tiny initial-condition errors grow with time and ultimately limit atmospheric predictability.

Model / Equation

Forecast Error Growth

The Lorenz model used in this section is identical to that introduced in Section 7.1. Throughout this section, the external forcing remains zero,

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

so that all model variability originates from the nonlinear Lorenz dynamics.

The climate model follows the same workflow introduced in Section 7.1.

Experiment configuration
        │
        ▼
Model initialization
        │
        ▼
Time integration
        │
        ├── Internal tendency (Lorenz model)
        ├── External forcing tendency
        ▼
Update model state
        │
        ├── Restart files
        ├── Model output
        ▼
Diagnostics and visualization

This workflow resembles the structure used by comprehensive climate models such as CESM and WRF. In this section, we reuse the same model framework and investigate how small differences in the initial condition influence the evolution of the simulated climate state.

To evaluate predictability, two simulations are performed.

The control simulation is denoted by

\[ \mathbf{X}_{ctl}(t) = \left( X_{ctl}, Y_{ctl}, Z_{ctl} \right), \]

while the perturbed simulation is denoted by

\[ \mathbf{X}_{prd}(t) = \left( X_{prd}, Y_{prd}, Z_{prd} \right). \]

Both simulations use

the same Lorenz model,
the same numerical integration,
the same model parameters, and
zero external forcing.

The only difference is the initial condition.

The perturbed simulation begins from

\[ \mathbf{X}_{prd}(0) = \mathbf{X}_{ctl}(0) + \delta\mathbf{X}_0, \]

where

\[ |\delta\mathbf{X}_0| \ll |\mathbf{X}_{ctl}(0)|. \]

Initially, the two simulations are nearly identical. As the integration proceeds, nonlinear dynamics amplify the small perturbation, causing the trajectories to separate.

The forecast error is defined as

\[ E(t) = \left| \mathbf{X}_{prd}(t) - \mathbf{X}_{ctl}(t) \right|. \]

The objective of this section is to investigate how forecast errors evolve with time and to compare their behavior between a simple linear model and the chaotic Lorenz model.

Sample Code

The sample programs in this section reuse the Model_src-v3.0 framework introduced in Section 7.1. 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_climate.py
│   ├── fig-timeseries_CTL_PRD-template.py
│   └── fig-error_growth-template.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 climate model.
`Lmodel_climate-v3.py`	Performs numerical integration of the Lorenz climate model.
`fig-timeseries_climate.py`	Generates a time-series plot for a single simulation.
`fig-timeseries_CTL_PRD-template.py`	Compares the control and perturbed simulations.
`fig-error_growth-template.py`	Calculates and visualizes forecast-error growth.

No new climate-model source code is introduced in this section. Instead, we reuse the same model framework developed in Section 7.1 and modify only the experiment setup and restart files to create control and perturbed simulations.

By comparing these two simulations, we investigate how very small initial-condition differences evolve in both a chaotic nonlinear system and a simple linear system.

Exercise

Step 1. Lorenz Control Simulation

First, perform the control simulation using the Lorenz model. This simulation will be used as the reference trajectory.

Create a new experiment directory,

mkdir Chapter8_1
cd Chapter8_1

Create a directory for the Lorenz control simulation,

mkdir CTL_Lorenz
cd CTL_Lorenz

Copy the model source code,

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

Open run_model_climate-v3.sh and set

REST="False"

TSTR="0"
TEND="100"

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

DEBUG_OUTPUT="False"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

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

ls out
ls rest

The restart file at model time 50 will be used as the initial condition for the perturbed simulation,

ncdump rest/restart.0050.nc

Next, create a common figure directory in the main experiment directory,

cd ..

mkdir figures

Copy the original time-series plotting script,

cp ../Model_src-v3.0/figures/fig-timeseries_climate.py \
   figures/fig-timeseries_CTL_Lorenz.py

Move to the figure directory,

cd figures

Open fig-timeseries_CTL_Lorenz.py and set the plotting period and the directory of control simulation output to

FDIR = "../CTL_Lorenz/out/"
TSTR = 40
TEND = 100

Run the plotting script,

python fig-timeseries_CTL_Lorenz.py

This figure shows the control trajectory from model time 40 to 100.

Questions

Why do we use the control simulation as the reference trajectory?
Why is the restart file at model time 50 useful for the next simulation?
What variables are shown in the control time-series figure?

Step 2. Lorenz Perturbed Simulation

Next, perform a second simulation that begins from a slightly perturbed initial condition. Rather than modifying the model equations or parameters, we introduce a tiny perturbation by editing one of the restart files generated during the control simulation.

Return to the main experiment directory,

cd ..

Create a directory for the perturbed simulation,

mkdir PRD_Lorenz
cd PRD_Lorenz

Copy the model source code,

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

Create a directory for the restart file,

mkdir rest

Copy the restart file corresponding to model time 50 from the control simulation,

cp ../CTL_Lorenz/rest/restart.0050.nc \
   ./rest/restart.0050.CTL.nc

Convert the restart file to CDL format,

ncdump rest/restart.0050.CTL.nc > restart.0050.cdl

Open restart.0050.cdl using a text editor and locate the model state variables,

time = 50 ;

Xint = 8.89485007017224 ;

Yint = 5.70490715205734 ;

Zint = 31.1349733070038 ;

Modify only the least significant digits of one or more variables. For example, round the restart values from

Xint = 8.89485007017224 ;
Yint = 5.70490715205734 ;
Zint = 31.1349733070038 ;

to

Xint = 8.89 ;
Yint = 5.70 ;
Zint = 31.13 ;

Generate the new restart file,

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

Verify that the restart file has been updated correctly,

ncdump rest/restart.0050.nc

Open run_model_climate-v3.sh and set

REST="True"

TSTR="50"
TEND="70"

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

DEBUG_OUTPUT="False"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

This simulation differs from the control simulation only through the tiny perturbation introduced into the restart file. The governing equations, model parameters, and external forcing are otherwise identical.

Next, make a simple time-series plot of the perturbed simulation.

Return to the main experiment directory,

cd ..

Copy the original plotting script into the common figure directory,

cp ../Model_src-v3.0/figures/fig-timeseries_climate.py \
   figures/fig-timeseries_PRD_Lorenz.py

Move to the figure directory,

cd figures

Open fig-timeseries_PRD_Lorenz.py and set

FDIR = "../PRD_Lorenz/out/"
TSTR = 50
TEND = 70

Run the plotting script,

python fig-timeseries_PRD_Lorenz.py

This figure shows the perturbed trajectory from model time 50 to 70.

Questions

What is the only difference between the control and perturbed simulations?
Why should the perturbation be much smaller than the magnitude of the model state?
Does the perturbed simulation look physically reasonable?
Before comparing with the control simulation, do you expect the two simulations to remain nearly identical throughout the integration? Explain your reasoning.

Step 3. Compare Simulations and Examine Error Growth

The perturbed simulation appears similar to the control simulation when viewed independently. To better understand the impact of the initial-condition perturbation, we first compare the two simulations on the same figure and then quantify their differences by calculating the error growth.

Return to the figures directory,

cd figures

Compare the Control and Perturbed Simulations

Copy the comparison template,

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

Open fig-timeseries_CTL_PRD.py and verify that the following settings are correct,

FDIR_CTL = "../CTL_Lorenz/out/"
FDIR_PRD = "../PRD_Lorenz/out/"

# Use a longer period for the control simulation
TSTR_CTL = 40
TEND_CTL = 80

# Use the prediction period for the perturbed simulation
TSTR_PRD = 50
TEND_PRD = 70

Run the plotting script,

python fig-timeseries_CTL_PRD.py

The resulting figure overlays the control and perturbed simulations. Notice that the two trajectories are nearly identical immediately after model time 50 but gradually separate as the integration proceeds.

Calculate and Visualize Error Growth

Next, quantify the differences between the two simulations.

Copy the error-growth template,

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

Open fig-error_growth.py and verify that the prediction period is

TSTR = 50
TEND = 70

Run the script,

python fig-error_growth.py

The figure displays the absolute error between the control and perturbed simulations,

|X_PRD − X_CTL|
|Y_PRD − Y_CTL|
|Z_PRD − Z_CTL|

using a linear y-axis.

Although the error clearly increases with time, it is difficult to determine whether the growth is exponential from a linear-scale figure. To better examine the growth rate, modify the script by changing the y-axis to a logarithmic scale.

For each subplot, add

for a in ax:
    a.set_yscale("log")

Run the script again,

python fig-error_growth.py

Compare the linear- and logarithmic-scale figures. On the logarithmic plot, the early stage of the error growth appears approximately as a straight line, indicating nearly exponential growth.

Questions

When do the control and perturbed simulations begin to diverge?
Which Lorenz variable exhibits the largest error growth?
Why is a logarithmic y-axis more suitable for visualizing error growth?
Does the error appear to grow approximately exponentially during the early stage of the prediction?
What does this result imply about the predictability of the Lorenz system?

Step 4. Linear Model Control Simulation

The Lorenz model demonstrated that a tiny initial-condition difference can grow rapidly because of nonlinear interactions among the model variables. To better understand the role of nonlinearity, we now replace the Lorenz equations with a simple linear model while keeping the same model framework.

Return to the main experiment directory,

cd ..

Create a directory for the linear control simulation,

mkdir CTL_Linear
cd CTL_Linear

Copy the Lorenz model source code,

cp ../CTL_Lorenz/Lmodel_climate-v3.py Lmodel_linear.py
cp ../CTL_Lorenz/run_model_climate-v3.sh .

Open Lmodel_linear.py.

First, replace

print(" == Run the Lorenz climate model ==")

with

print(" == Run the linear climate model ==")

Next, replace the Lorenz internal tendency,

def internal_tendency(X, Y, Z):

    dXdt = -SIGMA * X + SIGMA * Y
    dYdt = -X * Z + RHO * X - Y
    dZdt = X * Y - BETA * Z

    return dXdt, dYdt, dZdt

with

def internal_tendency(X, Y, Z, time):

    a = 5.0
    c = 2.5
    omega = 2.0 * np.pi / 12.0

    dXdt = -a * np.cos(omega * time)
    dYdt =  a * np.sin(omega * time)
    dZdt =  c * np.cos(omega * time)

    return dXdt, dYdt, dZdt

The new tendency depends only on time and no longer contains nonlinear interactions among the model variables.

Next, locate the model integration section and replace

dXdt, dYdt, dZdt = internal_tendency(Xint, Yint, Zint)

with

dXdt, dYdt, dZdt = internal_tendency(Xint, Yint, Zint, time[i])

Open run_model_climate-v3.sh and change the model name from

Lmodel_climate-v3.py

to

Lmodel_linear.py

Then set

REST="False"

TSTR="0"
TEND="70"

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

DEBUG_OUTPUT="False"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

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

ls out
ls rest

Next, create a time-series figure for the linear control simulation.

Return to the main experiment directory,

cd ..

Copy the plotting script,

cp ../Model_src-v3.0/figures/fig-timeseries_climate.py \
   figures/fig-timeseries_CTL_Linear.py

Move to the figure directory,

cd figures

Open fig-timeseries_CTL_Linear.py and set

FDIR = "../CTL_Linear/out/"
TSTR = 40
TEND = 100

Run the plotting script,

python fig-timeseries_CTL_Linear.py

The resulting figure shows a smooth periodic evolution of the linear model variables.

Questions

How does the linear control simulation differ from the Lorenz control simulation?
Which nonlinear terms were removed from the Lorenz equations?
Why do you expect this linear model to be more predictable than the Lorenz model?
Do the model variables exhibit periodic behavior?

Step 5. Linear Model Perturbed Simulation and Visualization

Finally, repeat the initial-condition sensitivity experiment using the linear model. The procedure is identical to that used for the Lorenz model, allowing a direct comparison between the nonlinear and linear systems.

Return to the main experiment directory,

cd ..

Create a directory for the perturbed linear simulation,

mkdir PRD_Linear
cd PRD_Linear

Copy the model source code from the control simulation,

cp ../CTL_Linear/Lmodel_linear.py .
cp ../CTL_Linear/run_model_climate-v3.sh .

Create a directory for the restart file,

mkdir rest

Copy the restart file corresponding to model time 50,

cp ../CTL_Linear/rest/restart.0050.nc \
   rest/restart.0050.CTL.nc

Convert the restart file to CDL format,

ncdump rest/restart.0050.CTL.nc > restart.0050.cdl

Open restart.0050.cdl using a text editor.

Locate the model state variables,

time = 50 ;

Xint = -8.28241453757486 ;

Yint = 14.7529867493495 ;

Zint = 4.14120726878743 ;

Modify only the least significant digits. For example, round the restart values to

Xint = -8.28 ;

Yint = 14.75 ;

Zint = 4.14 ;

Generate the new restart file,

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

Verify that the restart file has been updated correctly,

ncdump rest/restart.0050.nc

Open run_model_climate-v3.sh and set

REST="True"

TSTR="50"
TEND="70"

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

DEBUG_OUTPUT="False"

Run the model,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

Next, compare the control and perturbed simulations.

Return to the figures directory,

cd ../figures

Copy the comparison template,

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

Open fig-timeseries_CTL_PRD_Linear.py and modify the input directories,

FDIR_CTL = "../CTL_Linear/out/"
FDIR_PRD = "../PRD_Linear/out/"

TSTR_CTL = 40
TEND_CTL = 80

TSTR_PRD = 50
TEND_PRD = 70

Run the plotting script,

python fig-timeseries_CTL_PRD_Linear.py

Compare the control and perturbed simulations. Unlike the Lorenz model, the two trajectories remain nearly parallel throughout the prediction period.

Next, calculate the error growth.

Copy the error-growth template,

cp ../../Model_src-v3.0/figures/fig-error_growth-template.py \
   fig-error_growth_Linear.py

Open fig-error_growth_Linear.py and modify the input directories,

FDIR_CTL = "../CTL_Linear/out/"
FDIR_PRD = "../PRD_Linear/out/"

TSTR = 50
TEND = 70

Run the script,

python fig-error_growth_Linear.py

Because the linear-model errors remain very small and nearly constant, Matplotlib may automatically display the y-axis using offset notation (for example, +2.40 × 10⁻³). This notation is convenient for scientific plotting but can make it difficult to interpret the actual error values.

To display the full values on the y-axis, disable the automatic offset notation by adding

for a in ax:
    a.ticklabel_format(axis="y", style="plain", useOffset=False)

Run the script again,

python fig-error_growth_Linear.py

Although the offset notation has been removed, the error curves are still difficult to distinguish because the default y-axis range is much larger than the error variation.

To zoom in on the small changes, prescribe a narrow y-axis range by adding

ax[0].set_ylim(0.00240, 0.00242)
ax[1].set_ylim(0.00290, 0.00310)
ax[2].set_ylim(0.00120, 0.00121)

Run the script once more,

python fig-error_growth_Linear.py

The resulting figure clearly shows that the errors remain nearly constant throughout the prediction period. In contrast to the Lorenz model, the linear model exhibits little amplification of the initial-condition perturbation.

Questions

How does the perturbed linear simulation differ from the perturbed Lorenz simulation?
Do the control and perturbed linear simulations remain close throughout the prediction period?
Why is it useful to disable the automatic offset notation when plotting very small errors?
Why is it useful to adjust the y-axis range for this figure?
Compare the error-growth figures for the Lorenz and linear models. What is the most important difference?
What do these results suggest about the role of nonlinear dynamics in limiting predictability?

Key Points

The Lorenz model is highly sensitive to its initial conditions. Tiny perturbations introduced through the restart file grow rapidly because of the model's nonlinear dynamics.
A restart file allows a simulation to continue from any previous model state. By modifying only the least significant digits of the restart variables, we can investigate initial-condition sensitivity while keeping all other model settings unchanged.
Comparing the control and perturbed simulations reveals how initially similar trajectories diverge over time in a chaotic system.
The error-growth figure provides a quantitative measure of predictability by showing how differences between two simulations evolve with time.
Replacing the Lorenz equations with a simple linear model removes the nonlinear interactions responsible for chaotic behavior.
In the linear model, the initial perturbation remains nearly constant throughout the prediction period, demonstrating that linear systems are generally much more predictable than chaotic nonlinear systems.
This experiment illustrates one of the fundamental concepts of atmospheric science: even with a perfect numerical model, small uncertainties in the initial state can ultimately limit weather predictability.

Homework

Repeat the initial-condition sensitivity experiment using a different restart time.

Instead of perturbing the restart file at model time 50, use the restart file at model time 60. You may use the same control simulations generated during this exercise. Only the perturbed simulations need to be repeated.

Create a new perturbed Lorenz simulation starting from the restart file at model time 60.
Generate the Lorenz error-growth figure using fig-error_growth.py.
Create a new perturbed linear-model simulation starting from the restart file at model time 60.
Generate the linear-model error-growth figure using fig-error_growth_Linear.py.
Compare the Lorenz and linear-model error-growth figures and answer the following questions.

Questions

Does changing the restart time alter the overall error-growth behavior of the Lorenz model?
Does the linear model remain predictable when initialized from a different restart time?
Compare the Lorenz and linear-model error-growth figures. Which model exhibits rapid error amplification, and why?
What do these experiments suggest about the role of nonlinear dynamics in limiting weather predictability?

Submission

Submit the following:

fig-error_growth.png
fig-error_growth_Linear.png

along with brief written answers to the homework questions.