Skip to content

Chapter 10.2 Analysis Updates in Data Assimilation

Learning Objectives

After completing this chapter, students will be able to

  • Explain the purpose of data assimilation and its Bayesian interpretation.

  • Describe how the analysis combines the first guess and observations using an analysis weight.

  • Implement a simple analysis cycle using the Incremental Analysis Update (IAU) method.

  • Investigate how different analysis weights influence the analysis.

  • Evaluate analysis performance using root-mean-square error (RMSE).

  • Interpret the advantages of IAU over an instantaneous analysis update.

Concepts

Why Do We Need Data Assimilation?

In Chapter 9, we learned that observations contain measurement errors. Because the exact state of the atmosphere is unknown, we generated an ensemble of initial conditions by perturbing the observations. The resulting ensemble forecast represented the uncertainty in the future evolution of the atmosphere.

However, observations are not the only source of information about the current atmospheric state. Before each forecast cycle, numerical weather prediction systems also have a short-range forecast from the previous analysis, commonly called the first guess or background.

Therefore, at the beginning of a new forecast cycle, we have two independent estimates of the atmosphere:

  • Observation: direct measurement of the atmosphere, but containing observation errors.
  • First guess (background): a short-range model forecast, but containing forecast errors.

Neither estimate is perfect. Instead of choosing one over the other, data assimilation combines both sources of information to produce a better estimate of the atmospheric state.

Bayesian Interpretation

Data assimilation can be interpreted using Bayesian probability.

Suppose the observation estimates

20°C ± 2°C

while the model first guess estimates

18°C ± 3°C

Both estimates contain uncertainty, represented by their error distributions.

The objective of data assimilation is to estimate the most probable atmospheric state by combining the information from both the observation and the first guess.

This combined estimate is called the analysis.

Because both sources provide useful but imperfect information, the analysis is generally more accurate than either the observation or the first guess alone.

Forecast, Observation, and Analysis

The forecast cycle consists of three states:

Observation

The observation provides measurements of the true atmosphere.

Although observations contain valuable information, they also include random measurement errors.

First Guess (Background)

The first guess is a short-range forecast started from the previous analysis.

Because the numerical model is physically consistent, the first guess provides useful dynamical information about the atmosphere. However, model errors cause the forecast to deviate from the true state over time.

Analysis

The analysis combines the first guess and the observations to obtain the best estimate of the current atmospheric state.

The analysis then becomes the initial condition for the next forecast cycle.

Observation Error and Forecast Error

Both observations and model forecasts contain uncertainties.

  • Observation error represents uncertainty in the measurements.
  • Forecast error represents uncertainty in the model prediction.

Because these two sources of information are independent, combining them generally reduces the overall uncertainty. As a result, the analysis is expected to be closer to the true atmospheric state than either the observation or the first guess alone.

Incremental Analysis Update (IAU)

The conventional analysis applies the entire analysis increment instantaneously whenever a new observation becomes available. This can produce abrupt jumps in the model trajectory.

The Incremental Analysis Update (IAU) method distributes the same analysis increment gradually over several model time steps. This produces smoother transitions while maintaining the observational constraint.

Model and Equations

Perfect-Model Framework

In this chapter, we continue using the perfect-model framework introduced in the previous chapter.

The control simulation represents the true atmospheric state, which is known only because we are performing a numerical experiment. In real-world weather prediction, the true atmospheric state is unknown.

Synthetic observations are generated by adding random Gaussian errors to the true state,

\[ x_o=x_t+\varepsilon_o, \]

where

  • \(x_t\) is the true state,
  • \(x_o\) is the observation,
  • \(\varepsilon_o\) is the observation error.

The observation error is assumed to follow a Gaussian distribution,

\[ \varepsilon_o \sim N(0,\sigma_o^2), \]

where \(\sigma_o^2\) is the observation-error variance.

The first guess (background) is obtained by integrating the numerical model forward from the previous analysis,

\[ x_f=x_t+\varepsilon_f, \]

where

  • \(x_f\) is the first guess (forecast),
  • \(\varepsilon_f\) is the forecast error.

Similarly, the forecast error is assumed to be Gaussian,

\[ \varepsilon_f \sim N(0,\sigma_f^2), \]

where \(\sigma_f^2\) is the forecast-error variance.

Throughout this chapter, the observation and forecast errors are assumed to be unbiased and statistically independent.

Bayesian Estimation

At the beginning of a forecast cycle, we have two independent estimates of the atmospheric state:

  • the observation, containing observation errors, and
  • the first guess (forecast), containing forecast errors.

Neither estimate is perfect. Bayesian estimation combines these two sources of information to determine the most probable estimate of the true atmospheric state.

Using Bayes' theorem,

\[ P(x|y) = \frac{P(y|x)P(x)} {P(y)}, \]

where

  • \(P(x)\) is the prior probability (the first guess),
  • \(P(y|x)\) is the likelihood of the observations,
  • \(P(x|y)\) is the posterior probability.

The analysis corresponds to the state with the maximum posterior probability.

Assuming that the forecast and observation errors are unbiased, Gaussian, and statistically independent, the Bayesian estimate becomes a weighted average of the forecast and the observation.

Bayesian Analysis Equation

Let

  • \(T_f\) be the forecast (first guess),
  • \(T_o\) be the observation,
  • \(T_a\) be the analysis.

The Bayesian analysis is

\[ T_a = \frac{\sigma_f^2} {\sigma_f^2+\sigma_o^2} T_o + \frac{\sigma_o^2} {\sigma_f^2+\sigma_o^2} T_f, \]

where

  • \(\sigma_f^2\) is the forecast-error variance,
  • \(\sigma_o^2\) is the observation-error variance.

This equation has an intuitive interpretation.

If the observation error is much smaller than the forecast error,

\[ \sigma_o^2 \ll \sigma_f^2, \]

then the analysis is close to the observation.

Conversely, if the forecast error is much smaller than the observation error,

\[ \sigma_f^2 \ll \sigma_o^2, \]

then the analysis remains close to the forecast.

Thus, Bayesian estimation automatically assigns greater weight to the estimate with the smaller uncertainty.

Analysis Error

Combining two independent estimates reduces the uncertainty of the analysis.

The resulting analysis-error variance is

\[ \sigma_a^2 = \frac{\sigma_f^2\sigma_o^2} {\sigma_f^2+\sigma_o^2}. \]

Equivalently,

\[ \frac{1}{\sigma_a^2} = \frac{1}{\sigma_f^2} + \frac{1}{\sigma_o^2}. \]

Therefore,

\[ \sigma_a^2 < \sigma_f^2, \]

and

\[ \sigma_a^2 < \sigma_o^2. \]

The analysis uncertainty is smaller than both the forecast uncertainty and the observation uncertainty. Consequently, the analysis is expected to provide a more accurate estimate of the true atmospheric state than either estimate alone.

Analysis Equation for Data Assimilation

The Bayesian analysis can be expressed as a weighted average of the observation and the first guess,

\[ T_a = \frac{\sigma_f^2} {\sigma_f^2+\sigma_o^2}T_o + \frac{\sigma_o^2} {\sigma_f^2+\sigma_o^2}T_f. \]

Defining the analysis weight

\[ W = \frac{\sigma_f^2} {\sigma_f^2+\sigma_o^2}, \]

the Bayesian analysis equation becomes

\[ T_a = (1-W)T_f + WT_o. \]

Rearranging the equation gives

\[ T_a = T_f + W(T_o-T_f). \]

The term

\[ T_o-T_f \]

is called the innovation, representing the difference between the observation and the first guess.

The quantity

\[ W(T_o-T_f) \]

is called the analysis increment, which is added to the first guess to obtain the analysis.

Thus, the analysis equation can be interpreted as

\[ \boxed{ \text{Analysis} = \text{First Guess} + \text{Analysis Increment} } \]

This incremental form is widely used in operational data assimilation systems because it clearly separates the forecast from the observational correction.

In this chapter, the analysis weight is prescribed manually,

W = 1.0
W = 0.9
W = 0.5

to investigate how different levels of confidence in the observations influence the analysis.

In Chapter 10.2, we will estimate the forecast-error variance from an ensemble of forecasts and compute the analysis weight automatically using the Ensemble Kalman Filter (EnKF).

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 in the perfect-model framework. Synthetic observations are generated by adding random observation errors to the control simulation. A forecast initialized from the synthetic observations provides the first guess, and an analysis is obtained by combining the first guess with the observations using the Incremental Analysis Update (IAU) method.

Download the source-code archive,

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-perturbed.py
│   └── Lmodel_analysis-IAU-v3.py
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

The experiments use the following configuration:

Control simulation:      t = 0–100
Observation period:      t = 55–95
Observation error:       OERR = 1.0

The source files used in this chapter are summarized 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.
utils/get_observation.py Generates synthetic observations by adding random observation error to the control simulation.
utils/create_initial_obs-perturbed.py Creates an initial condition from a synthetic observation and optionally applies random perturbations to the initial state.
utils/Lmodel_analysis-IAU-v3.py Performs Lorenz model integration using the Incremental Analysis Update (IAU) method.
fig-timeseries_obs.py Plots the true state together with the synthetic observations.
fig-timeseries_CTL_PRD-template.py Provides a template for comparing the true state with a model forecast.

During the exercises, the original source programs are progressively extended to implement a simple data assimilation system based on the Incremental Analysis Update (IAU) method. Additional visualization and analysis programs are developed throughout the exercises and are therefore not included in the original Model_src-v3.0 source-code archive.

Exercise

Step 1. Generate the Control Simulation and Baseline 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 state for evaluating the analysis experiments developed later in this chapter.

Synthetic observations are then created by adding random observation errors to the control simulation. The baseline observation experiment uses

OERR = 1.0
OBS_DT = 1.0

The observation file is named

obs_err10_dt10.nc

where

err10    OERR = 1.0
dt10     OBS_DT = 1.0

This naming convention allows additional observation files with different error amplitudes and observation intervals to be stored in the same directory.

Create the Working Directories

Create the chapter and control-simulation directories,

mkdir Chapter10_1
cd Chapter10_1

mkdir CTL
cd CTL

Copy the model source files,

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=101
TSTEP=1

Run the control simulation,

bash run_model_climate-v3.sh

The model automatically creates the required output and restart directories and integrates the control simulation over

t = 0–100

Generate the Baseline Synthetic Observations

Copy the observation-generation script,

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

Open

get_observation.py

and configure

outfile_nc = "obs_err10_dt10.nc"

TSTR = 0
TEND = 100

OBS_TSTR = 45.0
OBS_TEND = 100.0

OERR = 1.0
OBS_DT = 1.0

Run

python get_observation.py

This creates

obs_err10_dt10.nc

The analysis experiments will begin at

t = 45

The period

t = 45–50

is treated as analysis spin-up, while the evaluation period is

t = 50–100

Visualize and Verify the CTL Simulation

Return to the chapter directory,

cd ../

mkdir Figure
cd Figure

Copy the figure script,

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

Modify

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

TSTR = 40
TEND = 100

Run

python fig-timeseries_obs.py

Verify that

  • the CTL trajectory is continuous,
  • observations begin at t = 45,
  • observations fluctuate around the true state, and
  • observations are separated by OBS_DT = 1.0.

Calculate the Variability of the True States

Open fig-timeseries_obs.py and remove the comment symbols (#) from the lines that read and plot the true states (Xtrue, Ytrue, and Ztrue). The true states will be used as the reference trajectory throughout this chapter.

Then add the following code after reading Xtrue, Ytrue, and Ztrue to calculate the standard deviation and variance of each variable.

print("========================================")
print("True state statistics")
print("========================================")

print( f"X: std = {np.std(Xtrue.values):8.3f}, " f"var = {np.var(Xtrue.values):8.3f}")
print( f"Y: std = {np.std(Ytrue.values):8.3f}, " f"var = {np.var(Ytrue.values):8.3f}")
print( f"Z: std = {np.std(Ztrue.values):8.3f}, " f"var = {np.var(Ztrue.values):8.3f}")

These statistics provide a useful reference for comparing the prescribed observation error with the natural variability of the Lorenz system.

Verify the Results

After completing this step, the directory should resemble

Chapter10_1/
├── CTL/
│   ├── out/
│   ├── rest/
│   ├── obs_err10_dt10.nc
│   ├── get_observation.py
│   ├── Lmodel_climate-v3.py
│   └── run_model_climate-v3.sh
│
└── Figure/
    └── fig-timeseries_obs.py

The CTL simulation represents the true evolution of the Lorenz climate system. The file obs_err10_dt10.nc provides the baseline synthetic observations used in the first analysis experiment. Later exercises will generate additional observation files with different observation errors and observation frequencies.

Step 2. Create and Run the Baseline Analysis with \(W=1.0\)

In this step, we develop the analysis model from the CTL model and perform the baseline analysis experiment using

OERR = 1.0
OBS_DT = 1.0
W = 1.0

The analysis equation is

\[ x_a=x_f+W(x_o-x_f). \]

When \(W=1.0\),

\[ x_a=x_o, \]

so the model state is replaced by the observation at each observation time.

The analysis begins at \(t=45\). The period from \(t=45\) to \(t=50\) is treated as analysis spin-up, and RMSE is evaluated over \(t=50\)–100.

Create the Baseline Analysis Experiment

Return to the chapter directory,

cd ../

Create the first analysis directory,

mkdir ANALYSIS10
cd ANALYSIS10

Copy the CTL model and run script,

cp ../CTL/Lmodel_climate-v3.py Lmodel_analysis-v2.py
cp ../CTL/run_model_climate-v3.sh run_model_analysis-v1.sh

Copy the initial-condition generator,

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

The directory should initially contain

ANALYSIS10/
├── create_initial_obs-perturbed.py
├── Lmodel_analysis-v2.py
└── run_model_analysis-v1.sh

Add the Analysis Settings to the Model

Open

Lmodel_analysis-v2.py

and add the following settings below the experiment settings:

# ============================================================
# Data-assimilation settings
# ============================================================
OBSFILE = "%OBSFILE"

WX = %WX
WY = %WY
WZ = %WZ

Read the Synthetic Observations

Before the time-integration section, add

# ============================================================
# Read observations
# ============================================================
if not os.path.exists(OBSFILE):
    raise FileNotFoundError(f"Observation file was not found: {OBSFILE}")

fin_obs = xr.open_dataset(OBSFILE)

time_obs = fin_obs.time
Xobs = fin_obs.Xobs
Yobs = fin_obs.Yobs
Zobs = fin_obs.Zobs

iobs = 0

while (iobs < time_obs.size and float(time_obs[iobs]) <= float(TSTR)):
    iobs += 1

The observation at \(t=45\) is already used as the initial condition, so the first analysis update uses the next available observation.

Add the Instantaneous Analysis Update

Locate the model-integration loop,

for i in range(NT):

Inside this loop, add the analysis update immediately before

# Save current total climate state

Use

# --------------------------------------------------------
# Analysis update
# --------------------------------------------------------
if iobs < time_obs.size:

    if np.isclose(
        float(time_obs[iobs]),
        float(time[i]),
        atol=DT / 2.0,
    ):

        # Current total forecast state
        Xf = Xint + Xext
        Yf = Yint + Yext
        Zf = Zint + Zext

        # Analysis increments
        dXa = WX * (float(Xobs[iobs].values) - float(Xf.values[0]))
        dYa = WY * (float(Yobs[iobs].values) - float(Yf.values[0]))
        dZa = WZ * (float(Zobs[iobs].values) - float(Zf.values[0]))

        # Apply the increments instantaneously
        Xint = Xint + dXa
        Yint = Yint + dYa
        Zint = Zint + dZa

        print( f"Analysis update: time={time[i]:.2f}" )

        iobs += 1

The analysis update must occur before the current state is stored in Xout, Yout, and Zout. Therefore, the output at an observation time represents the updated analysis state.

At the end of the model, close the observation file,

fin_obs.close()

Prepare the Initial Condition

Open

create_initial_obs-perturbed.py

and set

OBSFILE = "../CTL/obs_err10_dt10.nc"

Create the initial condition at \(t=45\) without an additional perturbation,

python create_initial_obs-perturbed.py --t_initial=45 --seed=1 --amp=0.0

This creates

rest/restart.0045.nc

Revise the Analysis Run Script

Open

run_model_analysis-v1.sh

and change the model source file from

MODEL_SRC="Lmodel_climate-v3.py"

to

MODEL_SRC="Lmodel_analysis-v2.py"

Configure the analysis period,

REST=True
TSTR=45
TEND=101
TSTEP=1

Set the observation file,

OBSFILE="../CTL/obs_err10_dt10.nc"

Set the baseline analysis weights,

WX=1.0
WY=1.0
WZ=1.0

Add the following substitutions to the sed command that creates model.py:

-e "s|%OBSFILE|${OBSFILE}|g" \
-e "s|%WX|${WX}|g" \
-e "s|%WY|${WY}|g" \
-e "s|%WZ|${WZ}|g" \

Run the baseline analysis,

bash run_model_analysis-v1.sh

The analysis is integrated from \(t=45\) through \(t=100\).

Create the Analysis Figure

Return to the Figure directory,

cd ../Figure

Copy the figure template,

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

Open

fig-timeseries_CTL_analysis1.py

and revise the input directories,

FDIR_CTL = "../CTL/out/"
FDIR_PRD = "../ANALYSIS10/out/"

Set the plotting and evaluation periods,

TSTR_CTL = 50
TEND_CTL = 100

TSTR_PRD = 50
TEND_PRD = 100

Revise the figure title,

figure_title = ("Control Simulation and Analysis (OERR = 1.0, W = 1.0)" )

Replace

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

with

ax[0].set_title(figure_title)

Change the prediction labels from

label="PRD"

to

label="Analysis"

Select the Common Evaluation Times

After reading the CTL and analysis output, select the CTL states at the exact analysis times,

Xtrue = Xctl.sel(time=time_prd)
Ytrue = Yctl.sel(time=time_prd)
Ztrue = Zctl.sel(time=time_prd)

Because both datasets use the same model time step, interpolation is not required.

Calculate RMSE

Add the following after reading the CTL and analysis data:

# ============================================================
# Calculate RMSE
# ============================================================
Xtrue = Xctl.sel(time=time_prd)
Ytrue = Yctl.sel(time=time_prd)
Ztrue = Zctl.sel(time=time_prd)

RMSE_X = np.sqrt(np.mean( (Xprd.values - Xtrue.values) ** 2))
RMSE_Y = np.sqrt(np.mean( (Yprd.values - Ytrue.values) ** 2))
RMSE_Z = np.sqrt(np.mean( (Zprd.values - Ztrue.values) ** 2))

print("========================================")
print("Analysis RMSE: t = 50–100")
print("========================================")
print(f"X RMSE = {RMSE_X:.3f}")
print(f"Y RMSE = {RMSE_Y:.3f}")
print(f"Z RMSE = {RMSE_Z:.3f}")

Add RMSE Values to the Figure

Add the RMSE value to each panel,

ax[0].annotate(
    f"RMSE = {RMSE_X:.3f}",
    xy=(0.12, 0.90),
    xycoords="axes fraction",
    horizontalalignment="right",
    backgroundcolor="white",
    fontsize=9,
)

ax[1].annotate(
    f"RMSE = {RMSE_Y:.3f}",
    xy=(0.12, 0.90),
    xycoords="axes fraction",
    horizontalalignment="right",
    backgroundcolor="white",
    fontsize=9,
)

ax[2].annotate(
    f"RMSE = {RMSE_Z:.3f}",
    xy=(0.12, 0.90),
    xycoords="axes fraction",
    horizontalalignment="right",
    backgroundcolor="white",
    fontsize=9,
)

Run the Figure Script

Run

python fig-timeseries_CTL_analysis1.py

This creates

fig-timeseries_CTL_analysis1.png

Verify that the figure shows

  • the CTL trajectory,
  • the baseline analysis trajectory,
  • the evaluation period \(t=50\)–100, and
  • RMSE values for \(X\), \(Y\), and \(Z\).

Verify the Results

After completing this step, the directory structure should resemble

Chapter10_1/
├── CTL/
│   ├── out/
│   ├── rest/
│   └── obs_err10_dt10.nc
│
├── ANALYSIS10/
│   ├── out/
│   ├── rest/
│   │   └── restart.0045.nc
│   ├── create_initial_obs-perturbed.py
│   ├── Lmodel_analysis-v2.py
│   └── run_model_analysis-v1.sh
│
└── Figure/
    ├── fig-timeseries_obs.py
    ├── fig-timeseries_CTL_analysis1.py
    └── fig-timeseries_CTL_analysis1.png

The ANALYSIS10 experiment provides the baseline analysis using

OERR = 1.0
OBS_DT = 1.0
W = 1.0

The resulting RMSE values will serve as the reference for the observation-error, analysis-weight, observation-frequency, and IAU experiments developed in the following steps.

Step 3. Sensitivity to Observation Error

In Step 2, we assumed nearly perfect observations using OERR = 1.0. In this step, we investigate how the quality of the analysis changes as the observation error increases while keeping the analysis weight fixed at W = 1.0.

Experiment Observation File OERR W
ANALYSIS10 obs_err10_dt10.nc 1.0 1.0
ANALYSIS10_ERR25 obs_err25_dt10.nc 2.5 1.0
ANALYSIS10_ERR50 obs_err50_dt10.nc 5.0 1.0

The observation interval remains OBS_DT = 1.0 for all experiments.

Create New Experiments

Return to the chapter directory,

cd ../

cp -r ANALYSIS10 ANALYSIS10_ERR25
cp -r ANALYSIS10 ANALYSIS10_ERR50

Remove the output directories copied from the baseline experiment,

rm -rf ANALYSIS10_ERR25/out
rm -rf ANALYSIS10_ERR50/out
Replace the copied restart files with the baseline restart file,

rm ANALYSIS10_ERR25/rest/*.*
rm ANALYSIS10_ERR50/rest/*.*

cp ANALYSIS10/rest/restart.0045.nc ANALYSIS10_ERR25/rest/
cp ANALYSIS10/rest/restart.0045.nc  ANALYSIS10_ERR50/rest/

All three experiments now begin from the same initial condition at t = 45. Therefore, any differences in the analysis results arise solely from the different observation errors.

Generate Additional Observation Files

Return to the CTL directory,

cd CTL

Open

get_observation.py

and change

OERR = 1.0

to

OERR = 2.5

Change the output observation file from

obs_err10_dt10.nc

to

obs_err25_dt10.nc

Run

python get_observation.py

to generate the new observation file.

Repeat the same procedure using

OERR = 5.0

and save the output as

obs_err50_dt10.nc

After completing this step, the CTL directory should contain

obs_err10_dt10.nc
obs_err25_dt10.nc
obs_err50_dt10.nc

These three observation files have the same observation interval (OBS_DT = 1.0) but different prescribed observation errors. They will be used to investigate the sensitivity of the analysis to observation uncertainty in the following experiments.

Configure and Run the Analysis Experiments

Move to the ANALYSIS10_ERR25 directory,

cd ../ANALYSIS10_ERR25

Open

run_model_analysis-v1.sh

and change the observation file from

obs_err10_dt10.nc

to

obs_err25_dt10.nc

Then run the analysis,

bash run_model_analysis-v1.sh

Next, move to the ANALYSIS10_ERR50 directory,

cd ../ANALYSIS10_ERR50

Again, open

run_model_analysis-v1.sh

and change the observation file to

obs_err50_dt10.nc

Finally, run the analysis,

bash run_model_analysis-v1.sh

After completing this step, the three analysis experiments

ANALYSIS10
ANALYSIS10_ERR25
ANALYSIS10_ERR50

will differ only in the prescribed observation error, allowing a direct comparison of the impact of observation uncertainty on the analysis.

Create the Figures

Copy the baseline plotting script:

cd ../Figure

cp fig-timeseries_CTL_analysis10.py fig-timeseries_CTL_analysis10_ERR25.py
cp fig-timeseries_CTL_analysis10.py fig-timeseries_CTL_analysis10_ERR50.py

Modify FDIR_PRD, the figure title, and the output filename for each experiment.

Run

python fig-timeseries_CTL_analysis10_ERR25.py
python fig-timeseries_CTL_analysis10_ERR50.py

Each script computes RMSE over t = 50–100.

Discussion

  1. How does RMSE change as observation error increases?
  2. Which variable is most sensitive to observation error?
  3. Why does larger observation error generally produce a poorer analysis?

Step 4. Sensitivity to the Analysis Weight

In the previous step, the analysis weight was fixed at

W = 1.0

which completely replaces the forecast with the observation at each analysis time.

In practice, observations contain errors, so it is often beneficial to retain part of the model forecast. In this step, we investigate how the analysis depends on the analysis weight while keeping the observation error fixed.

Experiment OERR W
ANALYSIS10_ERR25 2.5 1.0
ANALYSIS09_ERR25 2.5 0.9
ANALYSIS05_ERR25 2.5 0.5

The observation interval remains

OBS_DT = 1.0

for all experiments.

Create New Experiments

Return to the chapter directory,

cd ../

cp -r ANALYSIS10_ERR25 ANALYSIS09_ERR25
cp -r ANALYSIS10_ERR25 ANALYSIS05_ERR25

Remove the copied output files,

rm ANALYSIS09_ERR25/out/*.*
rm ANALYSIS05_ERR25/out/*.*

Replace the copied restart files with the baseline restart,

rm ANALYSIS09_ERR25/rest/*.*
rm ANALYSIS05_ERR25/rest/*.*

cp ANALYSIS10_ERR25/rest/restart.0045.nc ANALYSIS09_ERR25/rest/
cp ANALYSIS10_ERR25/rest/restart.0045.nc ANALYSIS05_ERR25/rest/

Configure and Run the Experiments

Move to the first experiment,

cd ANALYSIS09_ERR25

Open

run_model_analysis-v1.sh

and change

WX=1.0
WY=1.0
WZ=1.0

to

WX=0.9
WY=0.9
WZ=0.9

Run

bash run_model_analysis-v1.sh

Next move to

cd ../ANALYSIS05_ERR25

Again edit

run_model_analysis-v1.sh

and set

WX=0.5
WY=0.5
WZ=0.5

Run

bash run_model_analysis-v1.sh

Create the Figures

Return to the Figure directory,

cd ../Figure

Copy the baseline figure script,

cp fig-timeseries_CTL_analysis10_ERR25.py fig-timeseries_CTL_analysis09_ERR25.py
cp fig-timeseries_CTL_analysis10_ERR25.py fig-timeseries_CTL_analysis05_ERR25.py

For each figure script, modify

  • FDIR_PRD
  • the figure title
  • the output filename

Run

python fig-timeseries_CTL_analysis09_ERR25.py
python fig-timeseries_CTL_analysis05_ERR25.py

Each script computes RMSE over t = 50–100.

Discussion

Compare the three experiments.

  1. Does reducing the analysis weight improve or degrade the RMSE?
  2. Which weight produces the smallest RMSE?
  3. Why might a weight smaller than 1.0 perform better when observations contain errors?
  4. What happens if the analysis weight becomes too small?

Because the observations are imperfect, giving some weight to the model forecast may produce a more accurate analysis than completely replacing the forecast with every observation.

Step 5. Incremental Analysis Update (IAU)

In Step 4, the analysis-weight sensitivity experiments showed that using

W = 0.9

allows the analysis to follow the observations closely while retaining a small contribution from the model forecast.

However, the conventional analysis method applies the entire analysis increment at a single model time step. When the forecast and observation differ substantially, the analysis value can differ greatly from the model value one time step earlier, producing abrupt jumps or spikes in the analysis trajectory.

The Incremental Analysis Update (IAU) method reduces these abrupt transitions by distributing the analysis increment gradually over multiple model time steps.

In this exercise, we compare the conventional instantaneous analysis update with the IAU method while keeping the observation error, analysis weight, observation interval, initial condition, and model configuration unchanged.

Experiment OERR OBS_DT W Analysis method
ANALYSIS09_ERR25 2.5 1.0 0.9 Instantaneous update
ANALYSIS09_ERR25_IAU 2.5 1.0 0.9 Incremental Analysis Update

Only the analysis-update method changes between the two experiments.

Create the IAU Experiment

Return to the Chapter 10.1 directory.

cd ../

Confirm that the following directories are available:

ANALYSIS09_ERR25
Figure

Copy the conventional analysis experiment:

cp -r ANALYSIS09_ERR25 ANALYSIS09_ERR25_IAU

The new directory initially contains the same model code, shell script, restart file, and experiment settings as ANALYSIS09_ERR25.

Prepare the Directory

The copied directory contains the output and restart files from the conventional analysis experiment. Remove these files before running the IAU experiment.

rm ANALYSIS09_ERR25_IAU/out/*.*
rm ANALYSIS09_ERR25_IAU/rest/*.*

Next, copy the restart file used in the conventional experiment.

cp ANALYSIS09_ERR25/rest/restart.0045.nc  ANALYSIS09_ERR25_IAU/rest/

Both experiments now begin from the same model state at

t = 45

Therefore, any differences between the two experiments result solely from the analysis-update method.

Copy the IAU Model Code

The IAU model code is located in

Model_src-v3.0/utils/Lmodel_analysis-IAU-v3.py

Copy this file into the new experiment directory.

cp ../Model_src-v3.0/utils/Lmodel_analysis-IAU-v3.py ANALYSIS09_ERR25_IAU/

After copying, the directory should contain

Lmodel_analysis-v2.py
Lmodel_analysis-IAU-v3.py

The first file performs the conventional instantaneous analysis update, whereas the second file implements the Incremental Analysis Update (IAU).

Configure the IAU Run Script

Move to the experiment directory.

cd ANALYSIS09_ERR25_IAU

Open

run_model_analysis-v1.sh

Change

MODEL_SRC=Lmodel_analysis-v2.py

to

MODEL_SRC=Lmodel_analysis-IAU-v3.py

Do not change the observation file or the analysis weights.

The relevant portion of the shell script should now read

MODEL_SRC=Lmodel_analysis-IAU-v3.py

OBSFILE=../CTL/obs_err25_dt10.nc

WX=0.9
WY=0.9
WZ=0.9

Retain the same OBSFILE path used in ANALYSIS09_ERR25.

No other settings should be changed.

Run the IAU Experiment

Run the analysis model.

bash run_model_analysis-v1.sh

After the model finishes, confirm that new output files have been created.

ls out/

Create the IAU Time-Series Figure

Move to the Figure directory.

cd ../Figure

Copy the plotting script.

cp fig-timeseries_CTL_analysis09_ERR25.py \
   fig-timeseries_analysis09_ERR25_IAU.py

Open

fig-timeseries_analysis09_ERR25_IAU.py

Replace the conventional control directory

FDIR_CTL = "../CTL/out/"

with

FDIR_CTL = "../ANALYSIS09_ERR25/out/"

Change the analysis directory

FDIR_PRD = "../ANALYSIS09_ERR25/out/"

to

FDIR_PRD = "../ANALYSIS09_ERR25_IAU/out/"

Revise the legend labels so that the two trajectories are identified as

Instantaneous Analysis
IAU

Revise the figure title. For example,

plt.suptitle("Incremental Analysis Update (IAU): W = 0.9, OERR = 2.5")

Retain the plotting period

t = 50–100

and compute the RMSE using the same procedure as in the previous exercises.

Run the plotting script.

python fig-timeseries_analysis09_ERR25_IAU.py

The resulting figure directly compares the conventional instantaneous analysis with the IAU experiment. The conventional analysis may exhibit abrupt jumps at the observation times because the analysis increment is applied instantaneously. In contrast, the IAU trajectory should display smoother transitions because the same analysis increment is distributed gradually over multiple model time steps.

Key Points

  • The analysis combines information from both the first guess (forecast) and the observations using Bayesian estimation.

  • The analysis weight determines the relative influence of the first guess and the observations.

  • Greater confidence in the observations results in a larger analysis weight, while greater confidence in the first guess results in a smaller analysis weight.

  • The analysis uncertainty is generally smaller than either the forecast uncertainty or the observation uncertainty alone.

  • Using \(W=1\) completely replaces the first guess with the observation, whereas smaller values retain more information from the first guess.

  • Instantaneous analysis updates may introduce abrupt changes in the model trajectory.

  • Incremental Analysis Update (IAU) distributes the analysis increment over multiple model time steps, producing a smoother model evolution.

  • RMSE provides an objective measure of analysis accuracy and enables quantitative comparison among different analysis methods and analysis weights.

Homework

In this homework, investigate how the observation frequency influences the performance of the analysis.

Step 1. Change the Observation Frequency

In the exercises, observations are assimilated every model time unit using

OBS_DT = 1.0

Change the observation interval to

OBS_DT = 5.0

This reduces the observation frequency, allowing the model to evolve freely for a longer period between successive analysis updates.

Repeat the analysis experiments using the three prescribed analysis weights:

W = 1.0
W = 0.9
W = 0.5

Run the analysis cycle and generate the corresponding time-series and error figures for each experiment.

Step 2. Compare the Analysis Performance

Compare the results obtained with

OBS_DT = 5.0

with the original experiments using

OBS_DT = 1.0

Consider the following questions.

  • How does decreasing the observation frequency affect the analysis trajectory?
  • Does the analysis error increase when observations are less frequent?
  • Which analysis weight produces the smallest RMSE when OBS_DT = 5.0?
  • Does the best-performing analysis weight change when the observation interval increases?
  • Which variable (\(X\), \(Y\), or \(Z\)) is most sensitive to the reduced observation frequency?
  • Are the differences among \(W=1.0\), \(W=0.9\), and \(W=0.5\) larger when observations are less frequent?
  • Why does the model generally drift farther from the true state between less frequent analysis updates?

Submission

Submit

  • The time-series and error figures for the three experiments:
W = 1.0
W = 0.9
W = 0.5

using

OBS_DT = 5.0
  • A short paragraph (approximately 100–200 words) comparing the results with the original experiments using OBS_DT = 1.0.

Your discussion should explain how reducing the observation frequency affects the analysis accuracy and identify which analysis weight provides the best overall performance when observations are available every five model time units.