6.2 ESTIMATION OF MODEL ERRORS IN THE LOCAL ENSEMBLE TRANSFORM KALMAN FILTER
Hong Li1,*, Eugenia Kalnay1, Takemasa Miyoshi2 and Christopher M. Danforth3
1

Department of Atmospheric and Oceanic Science, University of Maryland, College Park, Maryland
2
Numerical Prediction Division, Japan Meteorological Agency
3
Department of Mathematics, University of Maryland, College Park, Maryland

1. INTRODUCTION1
Ensemble Kalman Filters (EnKF) have beenshown to be more accurate than 3D-Var in data
assimilation simulations under the assumption of aperfect model. However, in reality, the forecast
model has deficiencies and does not represent theatmospheric behavior precisely due to lack ofresolution, approximate parameterizations ofsubgrid scale physical processes, and numericaldispersion. For assimilation of real observations,
into an imperfect model, it is not yet clear whether
the EnKF will be competitive or better than thecurrent operational 3D-Var data assimilation
systems. Only recently have some EnKF schemes
advanced from the perfect model scenario to realworld situations. Houtekamer et al (2005) showedthat the quality of EnKF with perturbed
observations was comparable to 3D-Var. The
Ensemble Square-Root filter (EnSRF) was
reported (Whitaker et al 2004) to outperform the
NCEP 3D-Var in reconstructing the middle andlower tropospheric analysis in the Northern
Hemisphere at T62/L28 resolution. Miyoshi (2005)
also showed that EnKF is most advantageousover 3D-Var when the observing network is sparseand also that the advantage diminishes in thepresence of model errors. Model errors have astronger negative influence on the performance of
the EnKF than on the 3D-Var because Kalmanfiltering algorithms rely strongly on the assumptionof an unbiased model, an assumption which is notsatisfied in practice. Therefore, accounting forsystematic errors associated with model
deficiencies has become an important issue for alldata assimilation systems, and especially for
EnKF.
The Local Ensemble Transform Kalman Filter(LETKF) (Hunt 2005) is a relatively new data
assimilation scheme in the square root EnKF
family, and is similar to the Local Ensemble
Kalman Filter (Ott et al 2004) but faster. It has
been implemented to assimilate simulated
observations in the NCEP GFS model (Szunyogh
et al. 2005), and recently in the NASA fvGCM* Corresponding author address:

Hong Li 4339 CSS,
College Park, MD 20742; e-mail: Lhhelen@atmos.umd.edu

model (Liu et al. 2005). The results are excellent
for a perfect model scenario. In order to develop
the LETKF into a competitive, operationallyapplicable data assimilation system, it isnecessary to extend the application of the LETKF
to more realistic weather forecast systems byaccounting for systematic model errors.
The main goal in this work is to investigate
techniques for treating model errors in theEnsemble Kalman Filter, building on previous work,and to develop a data assimilation system capableof assimilating real weather observations. Though
we focus on the LETKF, the results may also be
applicable to other Ensemble Kalman Filter
methods. Recently, Baek et al (2005) extended
the work of Dee and Da Silva (1998), hereafterreferred as DDS, by using a high order bias
estimate scheme based on the state augmentation
correcting the model errors at each grid point, andalso accounting for the cross-correlation of model
state variables and bias. They successfully tested
this approach with the Lorenz 40-variable modeland showed the ability to correct forecast model
errors. Miyoshi (2005) tested the high order DDS
approach on the SPEEDY primitive equations
model (without cross correlation terms) with the
EnSRF assimilation scheme but found it to be
unsuccessful. He then tried a low order correction
approach developed by Danforth et al (2005) inwhich the model errors are expanded into loworder EOFs, with very good results. In this work,by using the LETKF as the data assimilation
method, we first test the Baek et al (2005) highorder approach with cross-correlations, todetermine whether the failure found by Miyoshi(2005) is due to the absence of cross-correlations.
However the high order approach of correcting at
each grid point is computationally very expensive.
Thus we develop a low order approach in the
LETKF, and compare both results, under the
hypothesis that model errors can be represented
by relatively few degrees of freedom and can thus
be efficiently corrected.

2. TECHNIQUESERRORSFOR

TREATING

MODEL

The major difference between (1), (2) and (5),(6) exists in the weighting matrix K b , the formerf

2.1 Full dimensional augmented state method

uses the bias forecast error covariance Pbb , theFriedland (1969) proposed the augmented
state method by augmenting the state vector by amodel bias vector to estimate both the state andbias variables. Building on this idea, Dee and Da
Silva (1998) developed a two-stage estimation
algorithm in which the estimation procedures forthe bias and the state are carried out successively.
The bias is estimated on every model grid point by
theobserved-minus-forecast
residuals

latter uses the cross covariance Pbx between thebias and state variables.

y o − h( x f )b a = b f − K b [ y o − h( x f ) + h(b f )](1)

where n is the time step, xn represents theanomaly of the state with respect to climatology,

K b = Pbbf H T (HPbbf H T + HPxxf H T + R) −1

(2)

f

f

f

where the matrix Pxx and Pbb are the forecast
error covariance for the bias and the stateis the operator
variables, respectively. h ·mapping the model state variables into
observation space.
The analysis for the state variables is thenobtained using the standard analysis procedure
with the unbiased forecast state ~x f = x f − ba .
Another augmented state method proposedby Baek et al (2005) is to obtain the optimal
estimates of the state and bias variables
simultaneously. They assume there are no bias
observations but that the bias on all model grids
can be updated from the state observations by thecross-correlation between the forecast state andbias through the cross covariance, while thiscross-correlation is assumed to be zero in DDS.
The unbiased state forecast is calculated firstfff~
as x = x − b , then is used in the augmented
state analysis procedure:⎛ xa ⎞ ⎛ ~xf⎞⎜ a ⎟ = ⎜ f ⎟ + K bx [ y o − h( ~x f )]⎜b ⎟ ⎜ b ⎟⎝ ⎠ ⎝ ⎠⎛P f
K bx = ⎜⎜ xxf⎝ Pbx

Pxbf ⎞⎛ H ⎞⎟⎜ ⎟⎟ ( HPxxf H T + R ) −1
f ⎟⎜Pbb ⎠⎝ 0 ⎠

(3)

T

(4)

Equations (3) and (4) implicate the bias is
estimated by

b a = b f + K b [ y o − h( ~x f )]= b f + K b [ y o − h( x f ) + h(b f )]K b = Pbxf H T ( HPxxf H T + R ) −1

2.2 Low order bias estimation scheme
Assume we have a simple nonlinear model such
that
m
xn+1
= Fnm (xnm ) = a m + b m xnm + c m xnm2 +… (7)

m

m

m

and a , b , c ,... are coefficients of stateindependent terms, and of terms that vary linearly,
quadratically or with higher order with the state
anomaly. These model parameters may be afunction of time or space, containing, for example,
a diurnal and an annual cycle. The model
approximates the evolution of a “true” system:

xnt +1 = Fnt (xnt ) = a t + b t xnt + c t xnt 2 + ... . (8)Here the super-indices m and t refer to model andtruth respectively. We define the state error (withthe sign of a correction) asδ xn = xnt − xnm ,so thatδ xn +1 = F (x ) − F (x )
t
n

t
n

m
n

m
n

(9)

We measure noisy observations of the truth

xnt = xnt + δ xnt
Then

δx n +1 = [δa + δbx nm + δcxnm 2 + ...](10)

+ [(b m + 2c m x nm 2 + ...)δx n ] + ο (δF )ο (δx n )tmtmHere δ F represents δa = a − a , δb = b − b ,δc = c t − c m ,… , the model parameter errors.
The first bracket on the RHS of (10) represents
error due to model deficiencies, i.e., the source of“external” error growth, and the second is thelinear tangent model or propagator of the error thatis the source of “internal” error growth related to
instabilities. We will neglect the last terms that areproportional to the product of parameter errors andstate errors.
If the true state was exactly known( δ xn +1 = 0 ), we could use the model (7) and initial
t

(5)

t

(6)

conditions xn and compute exactly the model
error at every time step. However, in practice weonly have noisy observations of the truth. If weassume that the observational error

δx

t
n

and

state error δ xn have zero mean and theirstandard deviation is small compared to the typicalsize of the model state anomaly, then we canmeasure equation (10) over many cases andobtainδ xn +1 = ⎡⎣δ a + δ bxnm + δ cxnm2 + ...⎤⎦ +randomerrorsIf we collect enough cases, we can obtain thesystematic errors due to model deficiencies bymlinear regression where the predictor is xn .
Danforth et al (2005) have performed thisprocedure with the SPEEDY model using 5 yearsof reanalysis as training for the linear regression.
Like Leith (1978) and DelSole and Hou (1999),
they retained only the state independent and thelinear terms. In the multiple linear regression, theyobtained the seasonally dependent model bias
and the diurnal errors that dominate the EOFs of
the state independent error corresponding to δ a .
The linear component corresponding to δ b was
obtained using SVD of the covariance of thecoupled model state anomalies and corresponding
measured errors with respect to the reanalysis.
This required spatial localization of the covariancematrix in order to reduce sampling problems. The
model bias δ a and the leading EOFs of the state
independent and state dependent errors explain
most of the model error and can be used as a loworder model correction.
In this study we plan to take advantage ofthese independently derived fields and assume
that during the data assimilation the actual errors
are proportional to them, with amplitudes that are
determined as part of the LETKF cycle.
(11)b =T βHere T denotes the pre-computed base fields
with dimension n×k, where each column denotes
base fields, including the seasonal bias, thediurnal errors corresponding to` δ a , and the state
dependent errors. The time-dependent amplitudesof these base fields are represented by thevector β (with dimension k) and can be estimatedin the analysis cycle. Since the dimension of β ismuch lower than the model dimension n,estimating β is far less costly than estimating b in
the analysis cycle.

3. IMPLEMENTATION ON THE LETKF IN THE
PRESENCE OF THE MODEL ERRORS
3.1 The SPEEDY model
The SPEEDY model (Molteni 2003) is a
recentlydeveloped
atmosphericgeneralcirculation model (AGCM) with simplified physicalparameterization
schemesthatarecomputationally efficient, but that maintain thebasic characteristics of a state-of-the-art AGCM
with complex physics. It has a spectral primitiveequation dynamics and triangular truncation T30
at 7 sigma levels.
3.2 LETKF data assimilation scheme
The Local Ensemble Transform Kalman Filter(LETKF) is chosen for this model error estimation
study since it belongs to the sequential data
assimilation family whose performance is sensitive
to model bias.
LETKF is an ensemble square-root filter inwhich the observations are assimilated to update
only the ensemble mean (shown in equation (12))
while the ensemble perturbations are updated bytransforming the forecast perturbations through atransform matrix (equation (13)) introduced by
Bishop et al (2001). The basic formulas used in
the LETKF (Hunt 2005) are given by~
x a = x b + X b P a ( HX b ) T R −1 [ y o − h( x b )] (12)~X a = X b [(k − 1) P a ]1 / 2
(13)
a

b

Here X , X
are the analysis and forecast
ensemble
perturbations,
respectively.
The~ a1/ 2transform matrix P

is the square-root of~a~amatrix (k − 1)P where P is given by−1~
P a = (k − 1) I + ( HX b ) T R −1 ( HX b )

[

]

(14)

with the dimension of k by k, where k is theensemble size, which is generally much smaller
than both the dimension of the model and thenumber of observations. Thus, the LETKF
performs the analysis in the space spanned by theforecast ensemble members, which greatly
reduces the computational cost. Furthermore,since the analysis is computed independently ateach grid point, the LETKF computation can be
performed in parallel.
3.3 Observations

The observations are obtained by adding thezero mean normal distributed noise to the
NCEP/NCAR reanalysis (NNR) fields (Kalnay et al1996), which are then an indirect estimate of thestate of the evolving atmosphere. With respect tothese observations, the SPEEDY model has
significant model errors. The observations areavailable on the model grid at every 4 grid points.
In each analysis cycle (6hour), the LETKF is
combined with one of the bias estimate methods
to correct the forecast errors due to inaccurate
initial conditions and model errors.
4. PRELIMINARY RESULTS

4.2 Effects of model errors on the LETKF
To investigate the effects of model errors, we
perform the LETKF using “realistic observations”
from NNR and “simulated observations” obtained
from a SPEEDY model “nature run” where nomodel errors exist. Figure 2 shows the strong
negative influence of the model errors on theperformance of the LETKF. In the presence ofmodel errors, the analysis is much worse than in
the perfect model scenario for u-wind, height, andother model variables (not shown) at all pressureslevels. Therefore, it is an important issue for
LETKF to estimate and correct the model errors.

4.1 Evidence of model bias
First we investigate the SPEEDY model bias
against the NNR. The SPEEY model is evolved
from the NNR every 6 hours. Figure 1 shows thedifferences between the SPEEDY 6hr forecasts
and the NNR verified at the same time, averaged
over two months in the period from January 1,
1982 to February 28, 1982, for the zonal wind andheight at 500 hPa. The largest model bias of the uwind can be seen in the polar regions. Orographic
effect is a major originator for the systematic
errors in the height field.

FIG.2 Global mean analysis RMS error at all pressurelevels averaged over the second month (Feb 1982) after
the initial transient behavior with LETKF in the perfectmodel scenario ( blue line) and in the presence of themodel errors ( red line) for u-wind (left) and height (right).

4.3 Constant mean bias correction
Following Miyoshi (2005), we simplify the loworder estimation scheme by correcting only thetime-averaged fields. In this way, T in equation(11) has just one column, given by the mean bias
over the two months period. β is fixed to 1.
Figure 3 illustrates the success of this simple bias
correction. Comparing the analysis RMSE at
500hPa with (green line) and without (red line) themean bias correction, we can see clearly that
subtracting the mean bias largely reduced theanalysis errors. This improvement is also
substantial at other pressure levels especially at
the lower levels (Figure 4).
FIG. 1 The SPEEDY 6-hour model bias of u-wind (unit:m/s) and height (unit: m) at 500hPa against NNR
reanalysis fields. The time averaged NNR over thesame period and at the same level is shown by contours.

FIG.3 Time evolution of the global mean analysis errorsof the u-wind (top) and height (bottom) fields at 500hPa
in the presence of the model errors ( red line) and in thecase of mean bias correction.

each grid point is computationally expensive andcan lead to very slow convergence of the LETKF
because it doubles or triples the number of
unknowns to be determined. According to Baek et
al (2005), in order to obtain the optimal analysis,
the number of the ensemble members also needs
to be doubled, making the computationally costeven higher. In contrast, the low order approachincreases the size of the model state by a verysmall percentage, the expectation is that it will bemuch more efficient than an augmentation with thefull model error field.
Our results have shown that simply
subtracting the constant mean bias from thebackground fields at every analysis cycle has asignificant positive impact on the LETKF. In our
future work, the time-variant amplitude of thismean bias field will be estimated by using the loworder method. We believe this will give us abetter result than the constant mean bias
subtraction. Moreover, we will correct the diurnal
bias and then the state-dependent bias. Retaining
each component of the bias estimate, we shouldbe able to get a good estimate of the true model
errors. If successful, this approach could not only
improve the performance of the LETKF, but also of
the forecast, as well as providing information onmodel errors useful for diagnostic purposes.

REFERENCES

FIG.4 Global mean analysis RMS at all pressure levels
averaged over the second month (Feb 1982) beyond
the spin-up time in the presence of the model errors( red line) and in the case of mean bias correction(green line) for u-wind (left) and height (right).

5. DISCUSSION AND FUTURE WORKThe full dimensional bias estimation method
has not been tested so far. We will test theaugmented state method in our future work.
However this high order approach of correcting at

Baek, S.J., B. R., Hunt, E. Kalnay, E. Ott, and I.
Szunyogh 2005: Local ensemble Kalmanfiltering in the presence of model bias,
submitted to Tellus.
Bishop, C. H., Etherton, B. J., and Majumdar, S. J.,
2001: Adaptive Sampling with the Ensemble
Transform Kalman Filter. Part I: Theoretical
Aspects.
Monthly Weather Review, 129,420-436.
Danforth, C., E. Kalnay, and T. Miyoshi, 2005:
Estimating and Correcting Global Weather
Model Error. In preparation.
Dee, D. P. and A. M. Da Silva, 1998: Dataassimilation in the presence of forecast
bias. Quart. J. Roy. Meteor. Soc., 124, 269295.
DelSole T. and A. Y. Hou, 1999: Empricalcorrection of a dynamical model. part 1:
Fundamental issues, Monthly. Weather.
Review, 127, 2533–2545.
Friedland, B. 1969: Treatment of Bias in Recursive
Filtering. IEEE Transaction on AutomaticControl, AC-14, 359-367.
Houtekamer, P.L, Mitchell, H.L., Pellerin, G.,
Buehner, M., Charron, M., Spacek, L., andHansen, B., 2005: Atmospheric data
assimilation with the ensemble Kalman filter:
Results with real observations. Monthly
Weather Review, 133, 604-620.
Hunt, B. R., 2005: Efficient Data Assimilation for
Spatiotemporal Chaos: a Local Ensemble
Transform Kalman Filter, In preparation.
Kalnay, E., M. Kanamitsu, R. Kistler, W. Collins, D.
Deaven, L. Gandin, M. Iredell, S. Saha, G.
White, J. Woollen, Y. Zhu, M. Chelliah, W.
Ebisuzaki, W. Higgins, J. Janowiak, K. C. Mo,
C. Ropelewski, J. Wang, A. Leetmaa, R.
Reynolds, Roy Jenne, and Dennis Joseph,1996: The NMC/NCAR 40-Year Reanalysis
Project. Bull. Amer. Meteorol. Soc, 77, 437–
471.
Leith, C. E., 1978: Objective methods for weather
prediction, Annual Review of Fluid Mechanics
10, 107–128.
Liu, J., E. Ferfig, H. Li, I. Szunyogh, B. R. Hunt, E.
Kalnay, E. Kostelich, and R. Todling, 2005:
Application of Local Ensemble Transform
Kalman Filter: perfect model experimentswith NASA fvGCM model. Submitted to the86th AMS Annual Meeting
Miyoshi, T. 2005: Ensemble Kalman filter
experiments with a primitive-equation global

model, PhD Thesis at the University of
Maryland.
Molteni, F., 2003: Atmospheric simulations using a
GCM with simplified physical parametrizations.
I: Model climatology and variability in multidecadal experiments. Climate Dyn., 20, 175191.
Ott, E., Hunt, B. R., Szunyogh, I., Zimin, A. V.,
Kostelich, E.J., Kostelich, M., Corazza, M.,
Sauer, T., Kalnay, E., Patil, D. J., and Yorke,
J. A., 2004: A Local Ensemble Kalman Filter
for Atmospheric Data Assimilation. Tellus,
56A, 415-428.
Szunyogh, I., E. J. Kostelich, G. Gyarmati, D. J.
Patil, B. R. Hunt, E. Kalnay, E. Ott, and J. A.
Yorke, 2005: Assessing a Local Ensemble
Kalman Filter: Perfect Model Experimentswith the NCEP Global Model. Tellus, 57A,
528-545.
Whitaker, J. S., G. P. Compo, X. Wei, and T.H.
Hammil,2004:
Reanalysis
without
Radiosondes
using
Ensemble
Data
Assimilation. Monthly Weather Review, 132,
1190-1200.