ABSTRACT
Title of dissertation: ENSEMBLE KALMAN FILTER EXPERIMENTS
WITH A PRIMITIVE-EQUATION GLOBAL MODEL
Takemasa Miyoshi, Doctor of Philosophy, 2005
Dissertation directed by: Professor Eugenia Kalnay
Department of Meteorology
The ultimate goal is to develop a path towards an operational ensemble Kalman
flltering (EnKF) system. Several approaches to EnKF for atmospheric systems have been
proposed but not systematically compared. The sensitivity of EnKF to the imperfections
of forecast models is unclear. This research explores two questions: 1. What are the
relative advantages and disadvantages of two promising EnKF methods? 2. How large
are the efiects of model errors on data assimilation, and can they be reduced by model
bias correction?
Chapter 2 contains a theoretical review, followed by the FORTRAN development
and testing of two EnKF methods: a serial ensemble square root fllter (serial EnSRF,
Whitaker and Hamill 2002) and a local EnKF (LEKF, Ott et al. 2002; 2004). We
reproduced the results obtained by Whitaker and Hamill (2002) and Ott et al. (2004) on
the Lorenz (1996) model. If we localize the LEKF error covariance, LEKF outperforms
serial EnSRF. We also introduce a method to objectively estimate the optimal covariance
in ation.
In Chapter 3 we apply the two EnKF methods and the three-dimensional varia-
tional method (3DVAR) to the SPEEDY primitive-equation global model (Molteni 2003),
a fast but relatively realistic model. Perfect model experiments show that EnKF greatly
outperforms 3DVAR. The 2-day forecast ?errors of the day? are very similar to the anal-
ysis errors, but they are not similar among difierent methods except in low ensemble
dimensional regions. Overall, serial EnSRF outperforms LEKF, but their difierence is
substantially reduced if we localize the LEKF error covariance or increase the ensemble
size. Since LEKF is much more e?cient than serial EnSRF when using parallel computers
and many observations, LEKF would be the only feasible choice in operations.
In Chapter 4 we remove the perfect model assumption using the NCEP/NCAR
reanalysis as the ?nature? run. The advantage of EnKF to 3DVAR is reduced. When we
apply the model bias estimation proposed by Dee and da Silva (1998), we flnd that the
full dimensional model bias estimation fails. However, if instead we assume that the bias
is low dimensional, we obtain a substantial improvement in the EnKF analysis.
ENSEMBLE KALMAN FILTER EXPERIMENTS WITH A
PRIMITIVE-EQUATION GLOBAL MODEL
by
Takemasa Miyoshi
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulflllment
of the requirements for the degree of
Doctor of Philosophy
2005
Advisory Commmittee:
Professor Eugenia Kalnay, Chair/Advisor
Professor James Carton
Professor Owen Thompson
Professor Brian Hunt
Dr. Joaquim Ballabrera
c Copyright by
Takemasa Miyoshi
2005
ACKNOWLEDGMENTS
First of all, I express my deep gratitude to Prof. Eugenia Kalnay for her advising
and supports during the entire processes of my graduate study. She always supported me
both academically and mentally, which was indispensable for me to complete this hard
task.
Secondly, I thank the Japanese Government, especially the National Personnel Au-
thority, the Ministry of Land, Infrastructure and Transport, and the Japan Meteorological
Agency, for supporting my graduate study through the Japanese Government long-term
fellowship program.
Thirdly, I truly appreciate Drs. Franco Molteni and Fred Kucharski for kindly
having provided the source codes and documents for the SPEEDY model. Their help was
essential in the present research.
Many thanks go to Dr. Joaquim Ballabrera and Christopher Danforth for their
useful comments and suggestions on the drafts of this thesis. I thank committee members:
Prof. Eugenia Kalnay, Prof. James Carton, Prof. Owen Thompson, Prof. Brian Hunt,
and Dr. Joaquim Ballabrera for their support.
I am grateful to Prof. Edward Ott and Dr. Istvan Szunyogh, for insightful dis-
cussions and kind support. John Harlim has developed LETKF based on the SPEEDY
model, corroborating my results of LEKF. I thank students in the chaos group?s weather
project, especially Dr. Shu-Chih Yang, John Harlim, Christopher Danforth, Junjie Liu,
and Hong Li, for many useful discussions. I thank Tomohide Higuchi for helpful discus-
ii
sions on some mathematical details of the recursive fllter theory. Kohei Aranami kindly
volunteered to let me use his powerful personal computer remotely. Some of the analy-
sis/forecast computations have been done on his computer. I am grateful to all my friends
who have supported my life in the U.S. in many ways.
Lastly, I thank my wife Seiko and my son Tetsushi for their support and patience
to encourage me even in living separately during my stay away from Japan.
Without tremendous help and support, I could not have achieved the hard task to
complete this thesis in the two-year graduate study. My sincere thanks for all the kindness
and supports are truly beyond description.
iii
TABLE OF CONTENTS
List of Tables viii
List of Figures xii
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 General methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 Theory of data assimilation methods and application on the Lorenz-96 model 9
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Theoretical review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 Optimal weighted mean . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Forecast error covariance . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Ensemble formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.4 Serial ensemble square root flltering (serial EnSRF) . . . . . . . . . 18
2.2.5 Local ensemble Kalman flltering (LEKF) . . . . . . . . . . . . . . . 22
2.2.6 Observational error covariance localization . . . . . . . . . . . . . . . 23
2.2.7 Covariance in ation . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.8 Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.3 Data assimilation experiments on a simple system . . . . . . . . . . . . . . 28
2.3.1 The Lorenz-96 model . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3.2 Implementation of data assimilation methods . . . . . . . . . . . . . 28
iv
2.3.3 Description of experiments . . . . . . . . . . . . . . . . . . . . . . . 29
2.3.4 Results in the case of 40 observations . . . . . . . . . . . . . . . . . 30
2.3.5 Results in the case of 20 observations . . . . . . . . . . . . . . . . . 34
2.3.6 Online estimation of covariance in ation parameter . . . . . . . . . . 42
2.3.7 Timing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.3.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3 Data assimilation on the SPEEDY primitive-equation model 46
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1.2 The SPEEDY model . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 3DVAR implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2.1 Theory for practical implementation . . . . . . . . . . . . . . . . . . 49
3.2.2 Implementation on the SPEEDY model . . . . . . . . . . . . . . . . 51
3.2.3 Background error statistics - NMC method . . . . . . . . . . . . . . 53
3.2.4 Response tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3 EnKF implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.1 Serial EnSRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.2 LEKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.4 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.1 Observational network . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.2 Data assimilation experiments . . . . . . . . . . . . . . . . . . . . . 67
3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5.1 3DVAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5.2 Serial EnSRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
v
3.5.3 LEKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.5.4 Timing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.5.5 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.6 Sensitivity experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.6.1 Response with a difierent observational network . . . . . . . . . . . 101
3.6.2 Usefulness of moisture observations . . . . . . . . . . . . . . . . . . . 103
3.6.3 Vertical error correlations . . . . . . . . . . . . . . . . . . . . . . . . 107
3.6.4 Error covariance localization . . . . . . . . . . . . . . . . . . . . . . 111
3.6.5 Random perturbation addition as ?stochastic seeding? . . . . . . . . 114
3.6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
3.7 Characteristics of the analysis and forecast errors . . . . . . . . . . . . . . . 116
3.7.1 Structures of the analysis increment . . . . . . . . . . . . . . . . . . 116
3.7.2 Analysis error and ensemble perturbation growth . . . . . . . . . . . 120
3.7.3 Error flelds and ensemble spread . . . . . . . . . . . . . . . . . . . . 123
3.7.4 Forecast errors, analysis errors, and bred vectors . . . . . . . . . . . 123
3.7.5 E-dimension and explained variance . . . . . . . . . . . . . . . . . . 128
3.7.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
3.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4 Data assimilation in the presence of model errors 138
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.1.2 OSSE methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.1.3 NCEP/NCAR reanalysis (NNR) . . . . . . . . . . . . . . . . . . . . 139
4.2 Model error statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
vi
4.2.1 Model error biases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.2.2 EOFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
4.3 Efiects of model errors on data assimilation . . . . . . . . . . . . . . . . . . 146
4.4 Model bias estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
4.4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
4.4.2 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 158
4.5 Summary and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
5 Summary and future directions 173
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
5.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
A The model-independent core modules of serial EnSRF and LEKF 179
A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
A.2 Serial EnSRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
A.3 LEKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
B Recursive fllter technique (Purser et al. 2003a) 183
B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
B.2 nth-order recursive fllter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
B.3 4th-order recursive fllter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
Bibliography 190
vii
LIST OF TABLES
2.1 The 1-year mean analysis RMSE of each method with optimal parameter
settings on the Lorenz-96 model when the number of observations is 40.
Observational error is given as 1.0, and the model natural variability with-
out data assimilation shows RMSE of 6.70. Ensemble size for the EnSRF
and LEKF is flxed at 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2 The same as Table 2.1 but in the case of 20 observations and 10 ensemble
members. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3 Analysis RMSE using LEKF with 8 ensemble members on Lorenz-96 model
when the number of observations is 20. Enhanced variance in ation is
applied, and the local patch parameter is flxed as l = 6. The other two
parameters, l2 and ?, are changed. . . . . . . . . . . . . . . . . . . . . . . . 38
2.4 Analysis RMSE using LEKF with 8 ensemble members on Lorenz-96 model
when the number of observations is 20. Online estimation of variance in a-
tion parameter is applied. l2 is flxed to l2 = 0. In the cases with observa-
tional error covariance localization, l is flxed to l = 15. . . . . . . . . . . . . 43
2.5 Computational time in seconds on a Linux PC with an Intel Celeron 2.7GHz
processor for a 425-day forecast-analysis cycle on the Lorenz-96 model. . . . 44
3.1 Vertical levels of the SPEEDY model outputs. Sigma levels are also used
as model levels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 Observational error standard deviation . . . . . . . . . . . . . . . . . . . . . 65
viii
3.3 Analysis RMSE of 500hPa height fleld in meters using 3DVAR on the
SPEEDY model when the sparse observational network is applied. The
RMSE is temporally averaged for a month (112 samples) after the ini-
tial one-month spin-up period. ?3DVAR(full)? denotes the case that full
spatial dependence of background error standard deviation is considered.
?3DVAR(zonal mean)? and ?3DVAR(hmean)? denote the cases that zonally
and horizontally uniform background error standard deviation is consid-
ered, respectively. ?NO ASSIM? denotes the case that no data assimilation
is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.4 Analysis RMSE of 500hPa height fleld in meters using serial EnSRF on
the SPEEDY model when the sparse observational network is applied. The
RMSE is temporally averaged for a month (112 samples) after the initial
one-month spin-up period. For comparison, the analysis RMSE of 3DVAR
in the same period is 31.14. The ensemble size (NBV) is chosen as 10,
20, or 30, and the horizontal length scale of the Schur product is chosen
as integers from 1 to 6. ?D? denotes fllter divergence. 4% multiplicative
covariance in ation is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.5 AnalysisRMSEof500hPaheightfleldinmetersusingLEKFontheSPEEDY
model when the sparse observational network is applied. The RMSE is
temporally averaged for a month (112 samples) after the initial one-month
spin-up period. For comparison, the analysis RMSE of 3DVAR in the same
period is 31.14. The ensemble size (NBV) is chosen as 10, 20, or 30, and
the local patch parameter l is chosen as integers from 1 to 6. ?D? denotes
fllter divergence. 4% multiplicative covariance in ation is applied. . . . . . 86
ix
3.6 AnalysisRMSEof500hPaheightfleldinmetersusingLEKFontheSPEEDY
model when the sparse observational network is applied and the parameter
l2 is changed. The ensemble size (NBV) is chosen to be 20 or 30, the lo-
cal patch parameter is flxed to l = 2. The RMSE is temporally averaged
for a month (112 samples) after the initial one-month spin-up period. For
comparison, the analysis RMSE of 3DVAR in the same period is 31.14. . . . 94
3.7 Timing (min:sec) of serial EnSRF with 10, 20, and 30 ensemble members
on the SPEEDY model when the sparse observational network is applied. . 96
3.8 Timing (min:sec) of LEKF with 10, 20, and 30 ensemble members on the
SPEEDYmodelwhenthesparseobservationalnetworkisapplied. Thelocal
patch parameter l is chosen as an optimal setting shown in the parenthesis,
and l2 is flxed to 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.9 Optimal localization parameters (length scale  of serial EnSRF and patch
size l of LEKF) with given ensemble sizes. Here, the covariance in ation
parameter is flxed to 4%, the parameter l2 in LEKF is flxed to 0. . . . . . . 97
3.10 Comparison of analysis RMSE of 500hPa height fleld in meters on the
SPEEDY model when the sparse observational network is applied. For
each ensemble size (NBV), best results are chosen for serial EnSRF and
LEKF from Tables 3.4 and 3.5. For 3DVAR, the best performance 31.14m
is chosen. ?NO ASSIM? denotes the case without data assimilation. . . . . 98
3.11 Timing (min:sec) of serial EnSRF and LEKF with 10, 20, and 30 ensemble
members on the SPEEDY model when the sparse observational network is
applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
x
3.12 Estimated timing (min:sec) of serial EnSRF and LEKF with 10, 20, and
30 ensemble members when 100 times more observations are assimilated
using a parallel computer with 100 computational nodes. Parallelization
e?ciency is assumed as 80% and 20% for LEKF and serial EnSRF, respec-
tively. We assume the increased observations require 50 and 100 times more
computation for LEKF and serial EnSRF, respectively. . . . . . . . . . . . 100
3.13 Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble mem-
bers with various localization scale of observational error covariance local-
ization. The local patch parameter is chosen to be l = 1,2. . . . . . . . . . 112
4.1 Analysis RMSE of 500hPa height fleld using EnKF with 10 ensemble mem-
bers and low order bias estimation using various numbers of bias basic
flelds. ?Uncorrected? denotes the case without bias correction. ?Mean
bias?, ?EOF1?, and ?EOF2? denotes the cases with low-order bias correc-
tion using the mean bias fleld, using the mean bias fleld and flrst EOF, and
using the mean bias fleld and flrst and second EOFs, respectively. . . . . . . 171
A.1 Arguments of the subroutine enkf anal0 that computes Kalman gain. . . . . 180
A.2 Arguments of the subroutine enkf serial that computes analysis ensemble
perturbations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
A.3 Arguments of the subroutine enkf schur that computes the factor of the
Schur product for localization. . . . . . . . . . . . . . . . . . . . . . . . . . 181
A.4 Arguments of the subroutine lekf core. ndiml, nobsl, and m denote the
dimension of the local patch, the number of observations in the local patch,
and ensemble size, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . 182
xi
LIST OF FIGURES
2.1 The Gaussian-like weighting function (dashed line) and Gaussian function
(solid line) when length scale  is set to 10. . . . . . . . . . . . . . . . . . . 21
2.2 Deflnition of the local regions around a center point (C) for LEKF in the
Lorenz-96 model. The larger region denoted by l is the local patch of
LEKF, the smaller region denoted by l2 is the subset of the local patch that
is averaged to obtain the global analysis. Since l denotes the length of one
side, the size of the local patch becomes 2l +1, similarly for l2. . . . . . . . 30
2.3 Parameter dependence of the analysis RMSE of serial EnSRF with 10 en-
semble members on the Lorenz-96 model when the number of observations
is 40, cf. Fig.3(b) of Whitaker and Hamill (2002). The horizontal and ver-
tical axes show the localization length scale  and the covariance in ation
paramter ?, respectively. The minimum error 0.20 is observed when  = 5
and ? = 0.02. ?FILTER DIVERGENCE? denotes the region with RMSE
of more than 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.4 The similar as Fig.2.3 but using LEKF with 10 ensemble members using
the enhanced variance in ation. The horizontal and vertical axes show the
local patch parameter l and the enhanced variance in ation paramter ?,
respectively. l2 is flxed to 1. The minimum error 0.20 is observed when
l = 5 and ? = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.5 The same as Fig.2.3 but in the case of 20 observations and serial EnSRF
with 10 ensemble members. The minimum error 0.33 is observed when
 = 6 and ? = 0.04. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
xii
2.6 The same as Fig.2.5 but serial EnSRF with 20 ensemble members. The
minimum error 0.30 is observed when  = 6 and ? = 0.02. . . . . . . . . . . 36
2.7 The same as Fig.2.5 but serial EnSRF with 8 ensemble members. The
minimum error 0.34 is observed when  = 5 and ? = 0.04. . . . . . . . . . . 37
2.8 The same as Fig.2.4 but LEKF with 8 ensemble members and enhanced
variance in ation. l2 is flxed to 0. The minimum error 0.35 is observed
when l = 5 and ? = 0.03. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.9 The same as Fig.2.8 but LEKF with 8 ensemble members and multiplicative
variance in ation. The way to in ate covariance is the same as that of serial
EnSRF. The minimum error 0.37 is observed when  = 5 and ? = 0.05. . . 40
2.10 The same as Fig.2.9 but LEKF with 8 ensemble members, multiplicative
variance in ation, and observational error covariance localization. The way
to in ate covariance is the same as that of serial EnSRF. The local patch
parameters are flxed to l = 15 and l2 = 0. The horizontal axis shows the
length scale  of observational error covariance localization. The minimum
error 0.33 is observed when  = 4 and ? = 0.03. . . . . . . . . . . . . . . . . 41
3.1 Zonal mean of the regression coe?cients r1 of u. The vertical and horizontal
axes show height and latitude respectively. The coe?cients re ect statistical
strength of the geostrophic relationship. The larger values, i.e. stronger
geostrophic balance, are observed in mid-latitudes. . . . . . . . . . . . . . . 54
3.2 Background error standard deviation of ps (top panel) and u at the 4th level
(bottom panel) that are estimated using the NMC method. The units for
ps and u are Pa and m/s, respectively. Strong spatial dependence especially
in zonal structures is observed, but their noisiness indicates sampling errors. 55
xiii
3.3 Background error horizontal correlations. The top panel shows vertical
structure of u in y-direction, where the vertical axis shows sigma levels. The
bottom panel shows the latitudinal structure of ps in x-direction, where the
vertical axis shows latitudes. The horizontal axis shows horizontal length
scale in grid spacing. The shades show correlation. . . . . . . . . . . . . . . 56
3.4 Vertical background error correlations of u-wind (solid lines with white
circle), v-wind (dashed lines with black circle), temperature (short dashed
lines with white square), and moisture (long-short dashed lines with black
square). The correlations are averaged horizontally. . . . . . . . . . . . . . . 58
3.5 Spatial response of the 3DVAR in the case of a single observation with the
observational increment of 1.0m/s at mid-latitude (top panel in m/s) and
high-latitude (bottom panel in m/s) at the 4th level ( = 0.51). . . . . . . . 60
3.6 Inter-variable response of the 3DVAR in the case of a single observation
with the observational increment of 1.0m/s at mid-latitude at the 4th level.
Analysis increment of T (top panel in K) and ps (bottom panel in Pa) are
shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.7 Analysis increment of 3DVAR (top panel) and the difierence between truth
and flrst guess (bottom panel). The flelds are u-wind flelds (in m/s) at the
4th level. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.8 Deflnition of the local regions around a center point for LEKF in the
SPEEDY model. The larger region denoted by l is the local patch of LEKF,
the smaller region denoted by l2 is the subset of the local patch that is av-
eraged to obtain the global analysis. Since l denotes a half of a side, the
area of the local patch becomes (2l +1)?(2l +1), similarly for l2. . . . . . 63
xiv
3.9 A regularly distributed dense observational network. 1056 stations (23 per-
cents of all grid points) are located at one grid point out of every four grid
points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.10 A regularly distributed sparse observational network. 264 stations (5.7
percents of all grid points) are located at one grid point out of every 16 grid
points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.11 An irregularly distributed realistic observational network. 415 stations (9.0
percents of all grid points) are located mostly over continents in the north-
ern hemisphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.12 Schematic of data assimilation experiments on the SPEEDY model. ?DA?
denotes the data assimilation system. The ?initial? fleld is the output by the
SPEEDY model which is difierent from the ?analysis? fleld because the high
frequency components are smoothed out by the spectral transformation of
the model. The difierence between ?initial? and ?truth? is deflned as the
analysis error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.13 Analysis RMSE of 500hPa height fleld in the case of 3DVAR with the
dense observational network (solid line), sparse observational network using
the same background error statistics based on the dense network (dashed
line), sparse observational network using background error standard devia-
tion multiplied by 2.0 (dotted line), and sparse observational network using
background error statistics based on the sparse network (dot-dashed line). . 71
xv
3.14 Analysis RMSE of 500hPa height fleld in the case of the sparse observa-
tional network with full spatial dependence of background error standard
deviation (solid line), with horizontally uniform background error standard
deviation (dashed line), and with zonally uniform background error stan-
dard deviation (dotted line). . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.15 Horizontal structure of analysis RMSE of height fleld in the case of 3DVAR
with horizontally uniform background error standard deviation and the
sparse observational network. Shades and contours show analysis RMSE
(in meters) and mean analysis flelds (in meters), respectively. Dots indi-
cate observational locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.16 Horizontal structure of analysis error bias (top panel) and standard devia-
tion (bottom panel) of height fleld (in meters) in the same case as Fig.3.15.
Shades show analysis error bias (top panel) and standard deviation (bot-
tom panel), contours show mean analysis flelds. Dots indicate observational
locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.17 Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of 3DVAR with hori-
zontally uniform background error standard deviation and the sparse ob-
servational network. Shades and contours show analysis RMSE and mean
analysis flelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
xvi
3.18 Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using serial EnSRF with 10 (solid line), 20 (dashed line), and 30
ensemble members (dotted line). The localization length scale of the Schur
product is chosen as optimal. For comparison, the RMSE of 3DVAR is
shown as the short dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.19 Horizontal structure of analysis RMSE of height fleld, similar to Fig.3.15,
in the case of serial EnSRF with 20 ensemble members and localization
scale of 2.0. Shades show analysis RMSE (in meters), contours show mean
analysis flelds (in meters). Dots indicate observational locations. . . . . . . 80
3.20 Horizontal structure of analysis error bias (top panel) and standard devia-
tion (bottom panel) of height fleld (in meters), similar to Fig.3.16, in the
case of serial EnSRF with 20 ensemble members and localization scale of
2.0. Shades show analysis error bias (top panel) and standard deviation
(bottom panel), contours show mean analysis flelds. Dots indicate observa-
tional locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.21 Analysis RMSE at the all pressure levels temporally averaged for one month
(112 samples) after the initial one-month spin-up period in the case of the
sparse observational network using serial EnSRF with 10 (solid line), 20
(dashed line), and 30 ensemble members (dotted line). The four panels (a),
(b), (c) and (d) correspond to u-wind, height, temperature (T) and speciflc
humidity (q) flelds, respectively. The observational error standard devia-
tions are shown as thin solid lines wherever applicable. The localization
length scale of the Schur product is chosen as optimal. For comparison, the
RMSE of 3DVAR is shown as the short dashed line. . . . . . . . . . . . . . 83
xvii
3.22 Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of serial EnSRF with 20
ensemble members and localization scale of 2.0. Shades and contours show
analysis RMSE and mean analysis flelds, respectively. . . . . . . . . . . . . 85
3.23 Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using LEKF with 10 (solid line), 20 (dashed line), and 30 ensemble
members (dotted line). The local patch parameter l is chosen as optimal.
For comparison, the RMSE of 3DVAR is shown by the short dashed line. . 87
3.24 Horizontal structure of analysis RMSE of height fleld, similar to Figs.3.15
and 3.19, in the case of LEKF with 20 ensemble members and local patch
parameter l = 3. Shades show analysis RMSE (in meters), contours show
mean analysis flelds (in meters). Dots indicate observational locations. . . . 89
3.25 Horizontal structure of analysis error bias (top panel) and standard devi-
ation (bottom panel) of height fleld (in meters), similar to Figs.3.16 and
3.20, in the case of LEKF with 20 ensemble members and local patch pa-
rameter l = 3. Shades show analysis error bias (top panel) and standard
deviation (bottom panel), contours show mean analysis flelds. Dots indicate
observational locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
xviii
3.26 Analysis RMSE at all pressure levels temporally averaged for one month
(112 samples) after the initial one-month spin-up period in the case of
the sparse observational network using LEKF with 20 (solid line) and 30
ensemble members (dashed line). The four panels (a), (b), (c) and (d)
correspond to u-wind, height, temperature (T) and speciflc humidity (q)
flelds, respectively. The observational error standard deviations are shown
as thin solid lines wherever applicable. The localization length scale of the
Schur product is chosen as optimal. For comparison, the RMSE of 3DVAR
is shown as the short dashed line. . . . . . . . . . . . . . . . . . . . . . . . . 92
3.27 Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of LEKF with 20 ensem-
ble members and local patch parameter l = 3. Shades and contours show
analysis RMSE and mean analysis flelds, respectively. . . . . . . . . . . . . 93
3.28 Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using 3DVAR (short dashed line), serial EnSRF (solid line) and
LEKF (dashed line). The ensemble size for serial EnSRF and LEKF is 30.
Parameters are optimal, that is, best results are shown for each method. . . 98
3.29 Analysis RMSE of 500hPa height fleld in the case of the realistic observa-
tional network using 3DVAR (short dashed line) and serial EnSRF (solid
line). The ensemble size for serial EnSRF is 30, the localization scale is 3.0. 102
3.30 Spatial distribution of analysis RMSE of 500hPa height flelds (in meters)
in the case of the realistic observational network using 3DVAR (top panel)
and serial EnSRF (bottom panel). The ensemble size for serial EnSRF is
30, and the localization scale is 3.0. Dots indicate the observational locations.104
xix
3.31 Zonal structure of analysis RMSE of u-wind fleld (top panel) and temper-
ature fleld (bottom panel) in the case of 3DVAR with the realistic obser-
vational network. Shades and contours show analysis RMSE and mean
analysis flelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.32 The same flgure as Fig.3.31 but in the case of serial EnSRF with 30 ensemble
members and localization scale of 3.0. . . . . . . . . . . . . . . . . . . . . . 106
3.33 Analysis RMSE of the u-wind fleld (a), hight fleld (b), temperature fleld
(c) and speciflc humidity fleld (d) in the cases of serial EnSRF with 10 en-
semble members with and without humidity observations shown by broken
and solid lines, respectively. The thin solid line shows observational error
standard deviation wherever applicable. . . . . . . . . . . . . . . . . . . . . 108
3.34 Impact of the vertical correlation in the analysis RMSE of 500hPa height
fleld using LEKF with 10 (black lines) and 20 (blue lines) ensemble mem-
bers. The black solid line and black short-dashed line show the cases with
and without the vertical correlation, respectively. The blue dotted line and
blue long-dashed line show the cases with and without the vertical correla-
tion, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
3.35 The same flgure as Fig.3.34 but for the full vertical correlation. . . . . . . . 111
3.36 Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble mem-
bers. Here, observational error covariance localization is applied. The solid,
long-dashed and short-dashed lines correspond to no localization, with lo-
calization scales  = 1.0, and  = 0.5 respectively. . . . . . . . . . . . . . . 113
xx
3.37 Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble mem-
bers. Here, observational error covariance localization is applied. Dashed
and solid lines show the cases with and without random perturbation ad-
dition, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
3.38 Background error fleld (shaded) and analysis increment fleld (contour) of
surface pressure (ps) at an arbitrary time in the case of 3DVAR (upper
panel) and serial EnSRF with 30 ensemble members and the localization
scale of 3.0 (bottom panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
3.39 Analysis zonal RMSE of 500hPa u-wind flelds in the cases of 3DVAR, serial
EnSRF with 30 ensemble members and the localization scale of 3, and
LEKF with 30 ensemble members, the local patch parameter l = 3 and
covariance localization scale  = 1.5. . . . . . . . . . . . . . . . . . . . . . . 119
3.40 The growth of the initial analysis errors and ensemble perturbations in the
cases of 3DVAR and EnKF. Thin solid line shows the growth of the initial
3DVAR analysis errors (RMSE(3DVAR)). Thick solid line and thick broken
line show the growth of the initial EnKF analysis errors (RMSE(EnSRF))
and ensemble perturbations (SPREAD(EnSRF)), respectively. . . . . . . . . 121
3.41 The error flelds (ananlysis/forecast minus truth; shaded) and the ensemble
spread of 500hPa height fleld (contour) at an arbitrary time. The top and
bottom panels show analysis and background (6-hour forecast), respectively,
both at the same time. Some parts have similar structures indicated by
rectangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
xxi
3.42 2-day forecast error fleld (shaded) and analysis error fleld (contour) of u-
wind at the 4th vertical level ( = 0.51) at an arbitrary time in the case
of serial EnSRF with 10 ensemble members and the localization scale of 1.
The correlation between the two flelds are 0.76. . . . . . . . . . . . . . . . . 125
3.43 The same as Fig.3.42, but in the case of serial EnSRF with 30 ensemble
members and the localization scale of 3. The correlation between the two
flelds are 0.79. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
3.44 The same as Fig.3.42, but in the case of 3DVAR. The correlation between
the two flelds are 0.74. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
3.45 Temporal mean of correlation coe?cients between 2-day forecast error and
analysis error flelds for four model variables (u, v, T, q) in the case of serial
EnSRF with 10 ensemble members and the local patch parameter l = 1. . . 127
3.46 5-day forecast error ((a) and (b)) and 2-day forecast error ((c) and (d)),
analysis error ((e) and (f)), E-dimension (g), and bred vector (h) of u-wind
at the 4th level all valid at the same time chosen arbitrarily, in the cases of
serial EnSRF with 10 members and 3DVAR. . . . . . . . . . . . . . . . . . 129
3.47 Zonal mean of the E-dimension (top panel) and explained variance (bottom
panel) in the case of serial EnSRF with 10 ensemble members temporally
averaged for one month (112 samples). . . . . . . . . . . . . . . . . . . . . . 131
3.48 Scatter plot of the zonal mean of the time-mean E-dimension and explained
variance, indicating high anti-correlation. . . . . . . . . . . . . . . . . . . . 132
3.49 The E-dimension (top panel) and explained variance (bottom panel) at the
4th level in the case of serial EnSRF with 10 ensemble members temporally
averaged for one month (112 samples). . . . . . . . . . . . . . . . . . . . . . 133
xxii
3.50 Scatter plot of the time mean E-dimension and explained variance at the
4th level, indicating high anti-correlation -0.66. . . . . . . . . . . . . . . . . 134
3.51 The E-dimension (top panel) and explained variance (bottom panel) at the
4th level in the case of serial EnSRF with 10 ensemble members at an
arbitrary time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
3.52 Scatter plot of E-dimension and explained variance at the 4th level at an
arbitrary time, indicating anti-correlation -0.26. . . . . . . . . . . . . . . . . 136
4.1 SPEEDY model biases of u, v [m/s] and T [K] relative to the NCEP/NCAR
reanalysis flelds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.2 SPEEDY model biases of q, z [m] and ps [Pa] relative to the NCEP/NCAR
reanalysis flelds. Each vertical level has difierent scales for q flelds, it is
omitted. Red color shows larger values. . . . . . . . . . . . . . . . . . . . . 143
4.3 EOFs of two-month model error samples from January 1, 1982 to February
28. The top and bottom panels show temperature flelds at the bottom level
(shades) and surface pressure flelds (contour) of the flrst and second EOFs,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.4 Time series of expansion coe?cients of the flrst to fourth EOFs for a week.
Solid, dashed, dotted, dash-dotted lines show flrst, second, third, and fourth
EOFs, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.5 Time series of expansion coe?cients of the flrst EOF (top panel) and the
third (solid line) and fourth (dotted line) EOFs (bottom panel) for 50 days. 149
xxiii
4.6 Analysis RMSE of 500hPa height fleld of 3DVAR (short dashed line) and
EnKF with 10 ensemble members (solid line) and 30 ensemble members
(long dashed line) in the case of the regular sparse observational network
(Fig.3.10) in the presence of model errors. . . . . . . . . . . . . . . . . . . . 150
4.7 Analysis RMSE of 500hPa height fleld of 3DVAR and EnKF with 30 ensem-
ble members in the case of the sparse and realistic observational networks in
the presence of model errors. Blue lines show the cases of realistic observa-
tional network, black lines show the cases of sparse observational network.
As for blue lines, short dashed line and long dashed line show 3DVAR and
EnKF, respectively. As for black lines, solid line and dotted line show
3DVAR and EnKF, respectively. . . . . . . . . . . . . . . . . . . . . . . . . 151
4.8 Analysis RMSE at the all pressure levels temporally averaged for one month
after the initial one-month spin-up period (112 samples) in the cases of
3DVAR (short-dashed lines), EnKF with 10 ensemble members (solid line),
EnKF with 30 ensemble members (long-dashed lines). Black and blue lines
show the cases with the regular sparse observational network and the real-
istic observational network, respectively. The four panels (a), (b), (c) and
(d) correspond to u-wind, height, temperature (T) and speciflc humidity (q)
flelds, respectively. The observational error standard deviations are shown
as thin solid lines wherever applicable. . . . . . . . . . . . . . . . . . . . . . 153
4.9 Analysis RMSE of 500hPa height fleld of EnKF with 10 ensemble members
in the case of the sparse observational networks in the presence of model
errors. The dashed and solid lines show the cases with and without Dee
and da Silva?s full dimensional bias estimation, respectively. . . . . . . . . . 159
xxiv
4.10 Analysis RMSE of 500hPa height fleld for EnKF with low order bias es-
timation (green short-dashed line), EnKF with constant bias subtraction
(black dotted line), 3DVAR with constant bias subtraction (blue dotted
line). Solid lines show the cases without bias estimation in the cases of
EnKF (black line) and 3DVAR (blue line), the same as solid and long-
dashed lines of Fig.4.6. Thin red line shows 6-hour forecast errors initiated
from the reanalysis flelds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.11 Values of all the seven components of the bias variable when the mean bias
fleld is used as a basis to span the bias where each vertical level is assumed
to be independent. The seven lines correspond to vertical levels. . . . . . . 162
4.12 Spatial distribution of analysis RMSE of 500hPa height flelds (shades in
meters) in the case of imperfect model experiments using EnKF with 10
ensemble members. Contours show mean analysis flelds (in meters). Top
and bottom panels show the cases without and with bias estimation, re-
spectively. Dots indicate the observational locations. . . . . . . . . . . . . . 164
4.13 Analysis error bias flelds of 500hPa height (shades in meters) in the case
of imperfect model experiments using EnKF with 10 ensemble members.
Contours show mean analysis flelds (in meters). Top and bottom panels
show the cases without and with bias estimation, respectively. Dots indicate
the observational locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
xxv
4.14 Analysis error standard deviation flelds of 500hPa height (in meters) in
the case of imperfect model experiments using EnKF with 10 ensemble
members. Contours show mean analysis flelds (in meters). Top and bottom
panels show the cases without and with bias estimation, respectively. Dots
indicate the observational locations. . . . . . . . . . . . . . . . . . . . . . . 166
4.15 Analysis RMSE at the all pressure levels temporally averaged for one month
after the initial one-month spin-up period (112 samples) in the cases of
EnKF with low order bias estimation (green short-dashed line), EnKF with
constant bias subtraction (black dotted line), 3DVAR with constant bias
subtraction (blue dotted line). Solid lines show the cases without bias
estimation in the cases of EnKF (black line) and 3DVAR (blue line), the
same as the black solid line and the black short-dashed line of Fig.4.8. The
four panels (a), (b), (c) and (d) correspond to u-wind, height, temperature
(T) and speciflc humidity (q) flelds, respectively. The observational error
standard deviations are shown as thin solid lines wherever applicable. . . . 168
4.16 Zonal structure of analysis RMSE of u-wind fleld (top panel in m/s) and
temperature fleld (bottom panel in K) in the case of imperfect model ex-
periments using EnKF with 10 ensemble members where no bias estimation
is applied. Shades and contours show analysis RMSE and mean analysis
flelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
4.17 The same flgure as Fig.4.16 but with bias estimation. . . . . . . . . . . . . . 170
xxvi
Chapter 1
Introduction
1.1 Background
In 2004, ten typhoons hit Japan, bringing tremendous damages to the local economy and
people?s lives. This was the worst typhoon season in 50 years. Better preparation for such
extreme weather would minimize damage and save lives. In order for people to prepare
well in advance of extreme events, it is important to predict typical weather events that
are directly connected to people?s lives as early and precisely as possible.
Modern weather forecasting is greatly dependent on numerical weather prediction
(NWP). By integrating momentum equations of the atmosphere with physical processes,
we can predict future state of the atmosphere. Since the problem is mathematically an ini-
tial value problem, forecasts require an estimate of the initial state of the atmosphere. In
the early stages of NWP, well-trained meteorologists created a subjective analysis by hand
to provide numerical values of the initial condition. Meteorologists would plot ?synoptic?
observations all over the world on a single map and draw a weather chart with constant
pressure lines, constant temperature lines, etc. However, even among well-trained me-
teorologists, the analyzed charts are difierent. For high quality NWP, it is essential to
generate an ?objective analysis? with less uncertainty given the same amount of informa-
tion.
Since the initial conditions are crucial to NWP, the objective analysis (data assim-
ilation) method has nearly as long a history as NWP itself, about 50 years. Many data
assimilation methods, such as successive corrections, optimal interpolation (OI), and the
1
three-dimensional variational method (3DVAR) have been used operationally as Daley
(1991) and Kalnay (2003) describe in their comprehensive textbooks. NWP models have
improved so much recently that they provide forecasts as precise as observations. As a
result, the efiect of errors in the initial condition is more important. The relative im-
pact of introducing a new data assimilation system has increased, and data assimilation
is recognized as one of the most important issues in the development of operational NWP
systems. Thus, more advanced data assimilation methods such as four-dimensional vari-
ational method (4DVAR) and Kalman flltering (KF) have been explored and some have
already been implemented operationally.
Modern methods of data assimilation are based on maximum likelihood estimation
using all available information, namely, observations and forecasts. Observations provide
direct measurements of the current state of the atmosphere, the model forecasts provide
information obtained from the past. Each source has limitations. Observations are subject
to errors in measurements as well as the limited number of observing stations. There is no
observational instrument without errors. Even if we have a perfect instrument, a grid point
value for an NWP model variable will be difierent from the perfect observation because
the grid point value represents a mean state for the surrounding area. The difierence is
known as the error of representativeness. In addition, we do not have observations for
all model grid points all over the world for all meteorological elements. Thus, we cannot
have perfect information for the present atmospheric state. As for the model forecasts, no
model that simulates Mother Nature can be perfect. In addition, the atmospheric system
is nonlinear, which introduces dynamically chaotic behavior. Even if we have a perfect
model, a chaotic system ?forgets? past information in a relatively short time, so only
limited information can be used from the past. In other words, if a perfect model were
2
known and the system were linear, all past observational information could be transmitted
toward the present, and we could reach better and better estimates of the truth with time
because of the accumulation of information. Thus, the efiects of imperfect observations,
an imperfect model, and the nonlinearity of the system introduce di?culties in estimating
the atmospheric state.
Furthermore, an important problem in atmospheric data assimilation lies in the
large number of degrees of freedom of NWP models. Even if we know a best method
for data assimilation in theory, the large number of degrees of freedom can prohibit its
explicit implementation in practice. Thus, a lot of simplifying assumptions need to be
made in atmospheric data assimilation. The simpliflcations will decrease as computational
capability increases.
Kalman flltering (KF, Kalman 1960) is a theory of data assimilation that optimally
combines all the available information for linear systems. The large number of degrees of
freedom of NWP models prohibit the direct implementation of KF. However, the recent
improvement of computational capability enables its implementation with simpliflcations.
In the ensemble formulation of Kalman flltering, a.k.a. ensemble Kalman flltering (EnKF)
initially proposed by Evensen (1994), a low rank of the forecast error covariance is assumed;
limited ensemble members are assumed to be able to represent forecast errors. After
Evensen (1994), several methods of EnKF have been proposed, the methods can be divided
into two groups: a perturbed observation (PO) method and a square root fllter (SRF)
method. The PO method uses an independent data assimilation cycle for each ensemble
member. Unfortunately, this method requires perturbing observational data (Burgers et al.
1998), a source of sampling errors. Alternatively, SRF does not introduce such sampling
errors, which is why Whitaker and Hamill (2002) suggested that the SRF is expected to
3
outperform the PO method.
Several SRF methods have been proposed: an ensemble transform Kalman fllter
(ETKF, Bishop et al. 2001), an ensemble adjustment Kalman fllter (EAKF, Anderson
2001), a serial EnSRF (Whitaker and Hamill 2002), all of which are efiective when obser-
vational data are assimilated serially, and a local EnKF (LEKF, Ott et al 2002; 2004),
where observations are assimilated simultaneously. All of these methods were proposed
and tested separately; there has been little research thus far comparing difierent methods
of SRF. In addition, as Tippett et al. (2003) mentioned, the formulation of SRFs is not
theoretically unique, and it is not clear if there is any optimal choice among difierent
implementations.
Lastly, an important problem with EnKF in practice lies in the imperfection of NWP
models. KF is optimal when the forecasting model is linear and perfect, both of which
are not the case in the real NWP. The problem of nonlinearity can be treated by taking a
short time interval for the data assimilation cycle if the observations are frequent enough so
that the tangent linear assumption is valid. However, advanced data assimilation methods
with fewer assumptions tend to be more sensitive to imperfection of models, since they
are more strictly optimal under the assumptions (e.g. Miyoshi 2004a). Thus, it is very
important to consider model errors in EnKF. Dee and da Silva (1998) discussed a way to
treat model biases in data assimilation. Assuming that model biases are another set of
prognostic variables and that observational increments, i.e. forecast minus observation,
are observations of model biases, we can construct an additional data assimilation system
for model biases. It is di?cult to derive a time evolving model for biases, and Dee and da
Silva (1998) used a ?persistent? model that is just the identity.
4
1.2 Objectives
The ultimate goal of the present research is to develop a path towards an operational
ensemble forecast-analysis system that is expected to produce substantially improved
weather/climate information. EnKF provides a way to achieve an ensemble forecast-
analysis system. However, as mentioned in the previous section, at least four ways to
implement EnKF on NWP models have been suggested, but the relative advantages and
disadvantages among them are not clear. In addition, EnKF has been tested most exten-
sively in the perfect model scenario, there is little research testing the sensitivities of EnKF
to imperfection of forecast models. Thus, to achieve the goal, it is relevant to answer two
questions:
1. What are the relative advantages of difierent implementations of EnKF?
2. What is the efiect of model errors on EnKF, and can we handle them?
Whitaker and Hamill (2002) suggested SRF is expected to outperform the PO
method, thus, only SRF is investigated in the present research. In addition, SRF can
be divided into two groups: serial treatment of observations and simultaneous treatment
of observations. Thus, two most e?cient methods of each group, a serial EnSRF (Whitaker
and Hamill 2002) and LEKF (Ott et al. 2002; 2004), are chosen. The present research
aims to investigate the following questions:
1. What are the relative advantages and disadvantages between a serial EnSRF and
the LEKF?
2. How large are the efiects of model errors on data assimilation, and how does a model
bias estimation work in a chosen EnKF method?
5
1.3 General methodology
To investigate the questions described in the previous section, data assimilation cycle ex-
periments on a primitive-equation global model are performed. The model used in the
present research is the SPEEDY model (Simplifled Parameterizations, primitivE-Equation
DYnamics, by Molteni 2003). The SPEEDY model has simplifled physical processes, the
resolution is T30L7 corresponding to grid size 96?48?7 and sigma vertical coordinate
( = p/ps, where p and ps are pressure and surface pressure, respectively). The prognostic
variables are horizontal wind velocity components (u, v), temperature (T), speciflc humid-
ity (q), and surface pressure (ps). Molteni (2003) has shown that the SPEEDY model has
characteristics similar to state-of-the-art models, though it has larger systematic errors.
The data assimilation experiments in this thesis are based on observational systems
simulation experiments (OSSEs), where a true ?nature run? is assumed to be known,
and observations are simulated by adding random noise to the nature run according to
observational errors. To investigate the flrst question, we flrst assume a perfect model.
Under the perfect model assumption, the nature run is created by the model integration
after a 1-year spin-up period. Three data assimilation systems (3DVAR, serial EnSRF,
and LEKF) are developed and their performances are compared. 3DVAR is representative
of classical methods, and it provides a control for the performance of data assimilation.
In the perfect model cases, sensitivity to the parameter values are investigated, and the
parameters are adjusted. Three types of observational locations are tested; a regularly
distributed dense network, a regularly distributed sparse network, and an irregularly dis-
tributed realistic rawinsonde network (few over oceans and many over land in NH). From
this experiment, one of the EnKFs (serial EnSRF) is chosen for the latter experiments.
To investigate the second question, we use the NCEP/NCAR reanalysis (Kalnay et
6
al. 1996; Kistler et al. 2001) to generate ?observations?. The NCEP/NCAR reanalysis
provides a best estimate of the real atmosphere using a recent operational forecast model
and a 3DVAR data assimilation system. The SPEEDY model is less realistic than the
1995 NCEP model used in the reanalysis. Thus, we can expect NCEP/NCAR reanalysis
flelds to be an approximation representative of the true evolution of the atmosphere for
the SPEEDY model. We perform the same experiments as in the perfect model case using
two data assimilation systems, 3DVAR and serial EnSRF. From the experiments using
the NCEP/NCAR reanalysis as ?truth?, the efiects of the model errors become clear. We
estimate model bias by including bias variables, assuming the observational increment is
an observation of the model bias and the forecast model of the model bias is persistence
(i.e. identity), following Dee and da Silva (1998).
1.4 Outline
Chapter 2 describes theory of data assimilation and experiments using the simple Lorenz-
96 model (Lorenz 1996; Lorenz and Emanuel 1998). Section 2.2 introduces the theory
of data assimilation in a self-contained manner. Data assimilation is introduced as an
optimal weighted mean between forecasts and observations. Then, we introduce forecast
covariance and the difierence between KF and OI. Ensemble formulation is introduced,
and serial EnSRF and LEKF are described. Before going to a realistic primitive-equation
atmospheric model, Section 2.3 describes data assimilation experiments under the perfect
model assumption using the simple Lorenz-96 model. This gives an idea of how each
method works and how it is difierent from each other. In addition, the results ensure the
program developed for the present research is correctly coded by comparing the perfor-
mance with those of Whitaker and Hamill (2002) and Ott et al. (2004). In this section, the
7
three methods (3DVAR/OI, serial EnKF and LEKF) and full KF are implemented and
tested on the Lorenz-96 model, the methods are compared to each other. Here, full KF
means that all flve KF equations, introduced in Section 2.2, are solved explicitly without
the low-rank approximation that is essential in EnKF.
Chapter 3 discusses the flrst question under the perfect model assumption. In this
chapter, 3DVAR, serial EnKF and LEKF are introduced separately, and implementations
of each method are described. After descriptions of implementation, the results are shown
and compared. Sensitivities to difierent parameters, observational networks, and local-
ization types are discussed. The characteristics of the analysis and forecast errors are
also discussed. For latter experiments, one EnKF method (serial EnSRF) is chosen in
this chapter. The choice is based on the performance (accuracy and e?ciency) that we
obtained using a single computer.
In Chapter 4, we remove the perfect model assumption. Using NCEP/NCAR reanal-
ysis flelds as the nature run, we introduce model errors. First, the model error statistics
are computed, accumulated by the difierence between the reanalysis flelds and the model
forecasts initiated by the reanalysis flelds. Then, data assimilation experiments assuming
the reanalysis flelds as truth are performed without model bias estimation, which clarifles
the efiects of the model errors on data assimilation. Finally, we apply the model bias
estimation scheme proposed by Dee and da Silva (1998) and discuss the results.
Chapter 5 summarizes the results and gives conclusion of the present research. This
chapter also describes possible future directions, based on our new perspective achieved
through the present research.
8
Chapter 2
Theory of data assimilation methods and application on the Lorenz-96 model
2.1 Introduction
We describe a theoretical review on data assimilation methods in Section 2.2. Section
2.3 describes data assimilation experiments on a simple system (the Lorenz-96 model,
Lorenz 1996; Lorenz and Emanuel 1998), providing an idea on how each method works
and how it is difierent from each other. Here, we develop the core modules used for this
entire research. The results shown in this chapter support the correctness of the Fortran
programs.
2.2 Theoretical review
This section describes a review on theory of data assimilation methods that are devel-
oped in the present research. Some basics in estimation theory are brie y introduced so
that the present dissertation is self-contained. For more comprehensive introduction and
details of atmospheric data assimilation, see the textbooks by Daley (1991) and Kalnay
(2003, Chapter 5). Bouttier and Courtier (1999) provide a more practical introduction
on atmospheric data assimilation. As for Kalman flltering, Jazwinski (1970) provides a
more mathematically precise introduction and discussion. Gelb et al. (1974) also provide
a more comprehensive and practical introduction on Kalman flltering. Bouttier (1996)
introduces Kalman flltering in the context of atmospheric data assimilation.
9
2.2.1 Optimal weighted mean
Data assimilation is a way to keep the model state close to the nature by assimilating
observations. Modern methods of data assimilation are a maximum likelihood estimation
using all available information, namely, forecasts xf and observations yo. We look for an
analysis state xa by a linear weighted mean between xf and yo:
xa = xf +K(yo ?H(xf)) (2.1)
where K denotes the weight (notation follows Ide et al. 1997). Since observations cannot
be taken at all grid points for all prognostic variables, the vectors yo and x lie in difierent
spaces. H converts prognostic variables x to the space of observational data yo and is
called an ?observational operator?. H is not necessarily a linear mapping.
Deflne difierences from the true state xt:
?xa = xa ?xt (2.2)
?xf = xf ?xt (2.3)
?yo = yo ?H(xt) (2.4)
as residual errors of analysis, forecast, and observation, respectively. Subtracting xt from
the both sides of eq.(2.1), we get
?xa ? ?xf +K(?yo ?H?xf) (2.5)
where H denotes the Jacobian of H. The tangent linear approximation made here is
accurate as long as ?xf is small. The analysis error covariance Pa is written as
Pa =
D
?xa(?xa)>
E
(2.6)
where h?i denotes the statistical expected value. Substituting eq.(2.5) into eq.(2.6), we
10
note
Pa =
??
?xf +K(?yo ?H?xf)
??
?xf +K(?yo ?H?xf)
?> 
=
D
(I?KH)?xf(?xf)>(I?KH)>
E
+
D
K?yo(?yo)>K>
E
= (I?KH)Pf(I?KH)> +KRK> (2.7)
where the forecast error covariance Pf and observational error covariance R are deflned
as
Pf =
D
?xf(?xf)>
E
(2.8)
R =
D
?yo(?yo)>
E
(2.9)
and the forecast error and observational error are assumed to be uncorrelated:
D
?xf(?yo)>
E
= 0 (2.10)
We assume the error probability distribution function (PDF) is Gaussian. Then the max-
imum likelihood state minimizes the total analysis error variance. That is, the optimal
weight gives the minimum of the trace of Pa. Thus,
?
?K (trace(P
a)) = 0 (2.11)
Using the formulas
?
?X
?
trace(XYX>)
?
= X(Y+Y>) (2.12)
?
?X (trace(XY)) = Y
> (2.13)
(cf. eqs.(2.1-72) and (2.1-73) of Gelb et al. 1974), we can transform eq.(2.11) with eq.(2.7)
as
?2(I?KH)PfH> +2KR = 0 (2.14)
11
Here, we used the fact that covariance matrices Pf and R are symmetric. By solving
eq.(2.14) for K, we get the optimal weight (a.k.a. Kalman gain)
K = PfH>(HPfH> +R)?1 (2.15)
For simplicity, deflne S as
S = HPfH> +R (2.16)
Substituting eq.(2.15) into eq.(2.7), we get the analysis error covariance
Pa = (I?PfH>S?1H)Pf(I?PfH>S?1H)> +PfH>S?1R(PfH>S?1)>
= Pf ?2PfH>S?1HPf +PfH>S?1HPfH>S?1HPf +PfH>S?1RS?1HPf
= Pf ?2PfH>S?1HPf +PfH>S?1(HPfH> +R)S?1HPf
= Pf ?2PfH>S?1HPf +PfH>S?1SS?1HPf
= Pf ?PfH>(HPfH> +R)?1HPf
= (I?KH)Pf (2.17)
Thus, the analysis error covariance is the forecast error covariance decreased by the factor
(I?KH).
2.2.2 Forecast error covariance
A forecast model M gives a mapping from an analysis state at time i (xai ) to a forecast
state at time i+1 (xfi+1)
xfi+1 = M(xai ) (2.18)
Deflne the mapping that gives the true time evolution to be Mt. Then
xti+1 = Mt(xti) (2.19)
12
A common assumption is that the model M contains random error
xti+1 = M(xti)?? (2.20)
where the random process ? has no bias but has covariance Q, that is,
h?i = 0 (2.21)
D
??>
E
= Q (2.22)
Then, the forecast error covariance matrix is written as
Pfi =
D
?xfi (?xfi )>
E
=
D
(xfi ?xti)(xfi ?xti)>
E
=
??
M(xai?1)?Mt(xti?1)
??
M(xai?1)?Mt(xti?1)
?> 
=
??
M(xai?1)?M(xti?1)+?
??
M(xai?1)?M(xti?1)+?
?> 
?
D?
M?xai?1 +???M?xai?1 +??>
E
=
D
M?xai?1(M?xai?1)> +M?xai?1?> +?(M?xai?1)> +??>
E
=
D
M?xai?1(?xai?1)>M>
E
+
D
??>
E
= MPai?1M> +Q (2.23)
where M (the Jacobian of M) is called the tangent linear model (TLM). The transpose
of TLM (M>) is called the adjoint model (ADJ). The tangent linear approximation is
valid as long as ?xa is small. In the derivation of eq.(2.23), the analysis error ?xa and the
random error ? are assumed to be uncorrelated, that is,
D
?xa?>
E
= 0 (2.24)
The flve equations (2.1), (2.15), (2.17), (2.18), and (2.23) constitute the Kalman
flltering algorithm. The optimal weight K is called the Kalman gain matrix in Kalman
13
flltering. In optimal interpolation (OI), eqs.(2.1) and (2.15) are used with Pf constant
in time. Usually, the time-independent forecast error covariance used in OI is denoted B
and is obtained from a long-term statistical mean.
2.2.3 Ensemble formulation
The covariance matrix P of size N ?N is deflned as
P =
D
?x(?x)>
E
(2.25)
? 1n?1
nX
l=1
?x(l)(?x(l))> (2.26)
where n and l denote the number of samples and a sample index respectively, and ?x is a
vector of length N (the number of degrees of freedom of the model). An inflnite number of
samples (n ! 1) is required to get a perfect estimate of P. Since the covariance matrix
P is real and symmetric, it has a real square root ?E:
P = ?E?E> (2.27)
Here, both P and ?E are N ?N matrices. The choice of ?E is not unique because ?EU also
satisfles eq.(2.27), where U is a unitary matrix, i.e., one that satisfles UU> = I. Eq.(2.27)
indicates that a given covariance can be perfectly represented by N samples
?E = 1p
N ?1
h
?x(1) ??? ?x(N)
i
(2.28)
Usually N is O(107) in NWP models, which prohibits the explicit estimation of the co-
variance matrix. However, it is common for the covariance matrix P to be degenerate (see
for example, Dee 1995; Fukumori and Malanotte-Rizzoli 1995; Cane et al. 1996). That
is, P has m << N non-zero eigenvalues. By ignoring zero eigenvalues, the eigenvalue
decomposition can be written as
14
P = SDS> (2.29)
= EE> (2.30)
where D is an m?m diagonal matrix whose components are m non-zero eigenvalues, and
S is an N ?m matrix whose columns are the corresponding eigenvectors. E is an N ?m
matrix, which can be written as
E = 1pm?1
h
?x(1) ??? ?x(m)
i
(2.31)
In this way, m samples express the covariance matrix P fairly well. As was the case for
?E, the choice of E is not unique because EU also satisfles eq.(2.30).
Using the square root of the covariance matrices, eq.(2.23) can be written as
Pfi = Efi (Efi )>
= MPai?1M>
= MEai?1(Eai?1)>M>
= MEai?1(MEai?1)> (2.32)
where the covariance of the model random process Q is ignored, i.e. we assume M = Mt.
Thus,
Efi = MEai?1 (2.33)
The TLM M gives the time evolution of the perturbation vector ?x. By the deflnition of
the TLM, the original nonlinear model can approximate eq.(2.33), thus,
MEai?1 = 1pm?1
h
M?xa(1)i?1 ??? M?xa(m)i?1
i
? 1pm?1
h
M(xai?1 +?xa(1)i?1)?M(xai?1) ??? M(xai?1 +?xa(m)i?1 )?M(xai?1)
i
? 1pm?1
h
M(xa(1)i?1)? ?xfi ??? M(xa(m)i?1 )? ?xfi
i
(2.34)
15
where xa(l) and ?xf denote the lth ensemble member and the ensemble mean, respectively:
xa(l)i = ?xai +?xa(l)i (2.35)
?xai = 1m
mX
l=1
xa(l)i (2.36)
?xfi = 1m
mX
l=1
M(xa(l)i?1) (2.37)
Eq.(2.34) is equivalent to the ensemble forecast with an ensemble of size m.
In the ensemble formulation, the optimal weight for data assimilation eq.(2.15) can
be written
K = PfH>(HPfH> +R)?1
= Ef(HEf)>[HEf(HEf)> +R]?1 (2.38)
Thus, we can avoid storing the entirePmatrix. We simply storeEf andHEf, which are m
members of the forecast ensemble and their equivalents in observational space respectively.
Since H> is not required, the nonlinear H can approximate the linearized H in a similar
way as in the nonlinear model (eq.(2.34)):
HEf = 1pm?1
h
H?xf(1) ??? H?xf(m)
i
? 1pm?1
h
H(xf +?xf(1))?H(xf) ??? H(xf +?xf(m))?H(xf)
i
? 1pm?1
h
H(xf(1))?H(?xf) ??? H(xf(m))?H(?xf)
i
(2.39)
Eq.(2.38) requires the inverse of a p ? p matrix where p denotes the number of
observations. If m < p, we can simplify computation of eq.(2.38) by expressing it as
K = Ef[I+(HEf)>R?1HEf]?1(HEf)>R?1 (2.40)
where the inverse of an m?m matrix is required. We assume here that difierent observa-
tions are uncorrelated to each other, therefore the observational error covariance matrix
R is almost diagonal, and its inverse is trivial.
16
Finally, the Kalman flltering algorithm requires the computation of Pa to obtain
Pf at the next analysis time. This process is equivalent to producing an appropriate
analysis ensemble in the ensemble formulation or ?ensemble update?. There are two types
of ensemble update methods:
1. a perturbed observation (PO) method
2. a square root fllter (SRF)
The PO method applies an independent data assimilation cycle to each ensemble member
yielding
Ea = (I?KH)Ef (2.41)
This gives
Pa = (I?KH)Pf(I?KH)> (2.42)
where the analysis covariance is too small compared to eqs.(2.7) and (2.17). Thus, this
process requires perturbing observational data so that we get the term KRK> in eq.(2.7)
(cf. Burgers et al. 1998). Whitaker and Hamill (2002) suggested that the PO method
has a disadvantage in perturbing observations because the perturbation introduces an
additional source of sampling errors.
Alternatively, SRF solves eq.(2.17) directly. Assume the ensemble update is given
by a linear combination of the forecast ensemble perturbations:
Ea = EfT (2.43)
Then, eq.(2.17) yields
Pa = EfTT>(Ef)>
= (I?KH)Ef(Ef)> (2.44)
17
By solving for T, we get an ensemble update formula. Since the choice of T is not unique
(TU also satisfles eq.(2.44)), the ensemble update formula is not unique. Among the
several methods of SRF, a serial method by Whitaker and Hamill (serial EnSRF, 2002)
and a batch method by Ott et al. (LEKF, 2002; 2004) are described in the present paper.
Tippett et al. (2003) and Evensen (2003) also review several methods of SRF.
In summary, the ensemble formulation of Kalman flltering with SRF is realized by
the following procedures:
1. Ensemble forecast provides Ef and ?xf (eq.(2.34))
2. Generate Kalman gain matrix K using eq.(2.38) or (2.40)
3. Compute the analysis ensemble mean ?xa using eq.(2.1) with xa and xf replaced by
?xa and ?xf
4. Get the analysis ensemble perturbation Ea by an ensemble update formula
5. Get the analysis ensemble members xa(l) using eq.(2.35)
6. The analysis ensemble members constitute initial conditions for the next ensemble
forecast; repeat the processes
2.2.4 Serial ensemble square root flltering (serial EnSRF)
The serial EnSRF (Whitaker and Hamill 2002) assumes an ensemble update of the form
Ea = (I? ?KH)Ef (2.45)
In order for eq.(2.45) to give a solution of eq.(2.17),
(I? ?KH)Pf(I? ?KH)> = (I?KH)Pf (2.46)
18
has to be satisfled (cf. eq.(9) in Whitaker and Hamill 2002). In this algorithm, data
assimilation is performed only once on the ensemble mean using eq.(2.1) with the optimal
weight eq.(2.38), eq.(2.45) is used for the ensemble update. The solution of ?K is given by
?K = PfH>h(HPfH> +R)?1/2i>h(HPfH> +R)1/2 +R1/2i?1 (2.47)
(Andrews 1968), cf. eq.(10) in Whitaker and Hamill (2002). This formula allows for
correlated observations. If observations are uncorrelated (R is diagonal) observations are
assimilated one at a time (serially), making the terms HPfH> and R scalar. In this case,
eq.(2.47) can be simplifled, and assuming ?K = fiK where fi is a scalar value, one can get
the formula for fi
fi =
?
1+
s
R
HPfH> +R
!?1
(2.48)
which was flrst derived by Potter in 1964, cf. (13) in Whitaker and Hamill (2002). Thus,
the explicit computation of ?K is avoided, and almost no additional computation is re-
quired. In the present research, observations are assumed to be independent of each
other, which makes only the computation of eq.(2.48) necessary. Thus, K and H are vec-
tors with N dimensions, R is a scalar. Observations are treated serially, hence the name
?serial EnSRF?.
In the ensemble formulation, we use only m ensemble members to reproduce the
N?N covariance matrixP, where m isO(102). This approximation is not usually justifled
because thePreproduced from the m members of ensemble contains large sampling errors,
resulting in artiflcially large covariances between distant points. Thus, it is necessary to
localize the in uence of observations. Intuitively, covariances between distant points can
be assumed to be zero, since it is reasonable to assume there is no correlation between, for
example, the temperature at College Park and the wind at Tokyo. In the serial EnSRF,
19
localization around observations is usually done using a Schur product with a weighting
function similar to Gaussian function but compactly supported (Gaspari and Cohn 1999;
Hamill et al. 2001). The Gaussian-like weighting function S ?f multiplies a factor S(r)
to an input function f (the correlation), where r denotes a measure of distance d from the
center point. In this way, the input function f is localized around the center point. The
factor S(r) is given by a flfth-order piecewise rational function:
S(r) =
8
>>>>
>>>
<
>>>
>>>>
:
1? 14r5 + 12r4 + 58r3 ? 53r2 (r ? 1)
1
12r
5 ? 1
2r
4 + 5
8r
3 + 5
3r
2 ?5r +4? 2
3r
?1 (1 < r ? 2)
0 (2 < r)
(2.49)
cf. eq.(14) in Hamill et al. (2001). Here, r is deflned as
r = dp10/3 (2.50)
where  denotes the characteristic length scale. Fig.2.1 shows how the weighting function
looks compared to the Gaussian function, given the same length scale  = 10. The
function has a shape similar to a Gaussian function, but it becomes exactly 0 for distant
points. In serial EnSRF, the Kalman gain K, now an N dimensional vector, is localized
around the observation point. In this way, we force the analysis increment to zero at
distant points from the observation, assuming the components of K in distant points are
caused by sampling errors. This method of ?gain localization? is equivalent to covariance
localization (cf. eq.(5) of Houtekamer and Mitchell 2001) as follows
K = (S ?Pf)H>
?
H(S ?Pf)H> +R
??1
(2.51)
? S ?
?
PfH>(HPfH> +R)?1
?
(2.52)
Here, we used the fact that HPfH is just a scalar in the serial treatment of observations.
If the observation is at a grid point, eq.(2.52) becomes exact.
20
?50 ?40 ?30 ?20 ?10 0 10 20 30 40 500
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1 Gaussian?like weighting function (sigma=10)
GaussianGaussian?like weighting
Figure 2.1: The Gaussian-like weighting function (dashed line) and Gaussian function
(solid line) when length scale  is set to 10.
21
2.2.5 Local ensemble Kalman flltering (LEKF)
LEKF (Ott et al. 2002; 2004; Szunyogh et al. 2005) treats the local patches surrounding
every grid point independently, which automatically avoids correlations between distant
points. Since the local region is much smaller than the globe, the degrees of freedom of
the local patch are much smaller than those of the globe, which is preferable for the low
rank assumption of the error covariance. Local patches corresponding to all grid points
are analyzed totally independently, and all the analyses are combined at the end to get
the global analysis. Hereafter, we consider only a local patch. m ensemble perturbations
span a space with r dimensions, where r ? m?1. The eigenvalue decomposition of Pf
Pf = G?PfG> (2.53)
deflnes a variable transformation G. ?Pf is a diagonal matrix with just r components
difierent from 0. Ignoring degenerate eigenvectors, we can deflne G and ?Pf as an N ?r
matrix and an r?r diagonal matrix, respectively. G> gives a transformation from the N
dimensional physical space to the r dimensional ?hat space?
?x = G>x (2.54)
?X = G>XG (2.55)
where x and X are an arbitrary N-dimensional vector and an arbitrary N ?N matrix,
respectively. The inverse transformation from the ?hat space? to the physical space is
given by
x = G?x (2.56)
In LEKF, physical vectors and matrices are transformed to the r-dimensional ?hat space?
where data assimilation is considered. After data assimilation processes, the analysis
22
ensemble in the ?hat space? is transformed back into the physical space. Analysis weighted
mean equations in the local patch are given by
?xa = ?xf + ?Pa ?H>R?1d (2.57)
?Pa = h(?Pf)?1 + ?H>R?1 ?Hi?1 (2.58)
= ?Pf
h
I+ ?H>R?1 ?H?Pf
i?1
(2.59)
(cf. eqs.(22) and (23) in Ott et al. 2004) which requires an inverse in the low dimensional
?hat space?. The ensemble update formula is given by
?Ea = ?EfY (2.60)
cf. eq.(34) in Ott et al. (2004), where Y is an m?m matrix that transforms the forecast
ensemble into the analysis ensemble in the ?hat space?. Eq.(2.60) yields
?Pa = ?EfYY>(?Ef)> (2.61)
cf. eq.(37) in Ott et al. (2004). As mentioned before, the choice of Y is not unique, but
an optimal choice of Y is given by
Y =
h
I+(?Ef)>(?Pf)?1(?Pa ? ?Pf)(?Pf)?1?Ef
i1/2
(2.62)
cf. eq.(41) in Ott et al. (2004). This choice of Y yields the minimum difierence between
the forecast ensemble perturbations and the analysis ensemble perturbations. This Y is
advantageous to the smoothness and physical balance of the analysis flelds.
2.2.6 Observational error covariance localization
In the LEKF formulation described in the previous section, LEKF localizes around an
analyzed point with Heaviside-like step weighting whereas serial EnSRF used Gaussian-
like weighting function. Hunt (pers. comm.) suggested a Gaussian-type error covariance
23
localization in LEKF. Originally, all observations had the same weight in the local patch.
However, when the ensemble size is very small, we expect sampling errors cannot be
negligible even with the small local patch size, and thus, farther away observations may
add large sampling errors. If observations are sparse compared to the local patch size,
there is an unfavorable situation that only a single observation far from a local patch
center is included in a patch. This may introduce a large sampling error in the analysis
at the local patch center.
Thus, Hunt suggested to approximate the localization of Pf (which is di?cult to
carry out in the LEKF formulation) by weighting observational error according to the
distance from the center of the local patch. Intuitively,
K = (??Pf)H>
?
H(??Pf)H> +R
??1
(2.63)
? PfH>
?
HPfH> +(??1 ?R)
??1
(2.64)
Here, ?? denotes the localization function applied to a state vector or a matrix whose
components have correspondence to grid points. ?? multiplies weights on the components
of the vector or matrix according to the distance for the localization. Thus, if we want
to localize the error covariance around the center of the local patch, we just multiply
inverse of the covariance localization weight to the observational error variance. If we
apply Gaussian-type weighting, we substitute the observational error variance R with
R ? R?exp
?
d2
2 2
!
(2.65)
where d and  denote distance from the center of the local patch and localization scale, re-
spectively. We call this localization ?observational error covariance localization? originally
suggested by Hunt (pers. comm.).
24
2.2.7 Covariance in ation
Kalman flltering applied in nonlinear systems usually underestimates the forecast error
covariance because of the nonlinearity of the system and model errors that we neglected
(eq.(2.33)). The forecast error covariance needs to be in ated for stable flltering, so that
observations are given enough weight. If not, the process amplifles feedback and eventually
results in fllter divergence.
A simple way to in ate covariance is to multiply the covariance matrix Pf or Pa by
a factor slightly larger than 1
Pf ?Pf ?(1+?) (2.66)
where ? is a small positive number (? << 1). Using the in ated covariance Pf, ob-
servations have larger efiects and we prevent fllter divergence. This operation is called
?multiplicative covariance in ation?.
Ott et al. (2004) proposed a more sophisticated way to in ate covariance in the
context of LEKF, which is called ?enhanced variance in ation?. The diagonal components
of the analysis error covariance matrix in the ?hat space? ?Pa are in ated as follows:
?Pa ? ?Pa + ? ?trace(?Pa)
r I (2.67)
where ? and r denote a small positive number and the rank of Pa, respectively (cf. eq.(42)
in Ott et al. 2004). In this way, the in ation magnitude depends on the total analysis
error variance. In addition, the in ation is the same for perturbation directions with any
uncertainty; an eigendirection corresponding to a small eigenvalue has the same in ation
as an eigendirection corresponding to a large eigenvalue.
25
2.2.8 Variational formulation
Assuming a Gaussian shape of the error statistics, we can get the probability density
function (PDF) with respect to the forecast state as
Probf(x) / exp
 
?12(x?xf)>B?1(x?xf)
?
(2.68)
In the same way, the PDF with respect to the observation is given by
Probo(x) / exp
 
?12(Hx?yo)>R?1(Hx?yo)
?
(2.69)
The joint probability is written as
Probb&o(x) = Probb(x)?Probo(x)
/ exp
 
?12(x?xb)>B?1(x?xb)? 12(Hx?yo)>R?1(Hx?yo)
?
(2.70)
Since data assimilation is a maximum likelihood estimate, the analysis state x is obtained
by maximizing eq.(2.70). This is equivalent to minimizing the cost function
J(x) = 12(x?xb)>B?1(x?xb)+ 12(Hx?yo)>R?1(Hx?yo) (2.71)
At a minimizer x?, the gradient of the cost function is equal to 0 , that is,
rJ(x?) = 0 (2.72)
3DVAR algorithms look for the solution that satisfles eq.(2.72) using a quasi-Newton
algorithm. The solution found by the quasi-Newton minimizer is not necessarily the
global minimum of the cost function. However, since we start from a state very close
to the global minimum, we assume the local minimum found coincides with the global
minimum. This assumption is reasonable considering we have a 6-hour model forecast
26
that is close to observations, and the most likely state lies in somewhere between the
forecast and observations.
In the incremental form, the cost function (eq.(2.71)) can be rewritten as
J(xb +?x) ? 12?x>B?1?x+ 12(H?x?d)>R?1(H?x?d) (2.73)
where we deflned analysis increment ?x and observational increment d as
?x = x?xb (2.74)
d = yo ?Hxb (2.75)
Eq.(2.73) is not an exact equation because of the tangent linear approximation of the
generally nonlinear operator H, which is valid as long as ?x is small. If the observational
operator is linear, eq.(2.73) becomes exact. The gradient of eq.(2.73) is
rJ(?x) = B?1?x+H>R?1(H?x?d) = 0 (2.76)
Solving eq.(2.76), we get
?x = (B?1 +H>R?1H)?1H>R?1d (2.77)
= BH>(R+HBH>)?1d (2.78)
which is equivalent to the weighted mean (eq.(2.1)) with the optimal weight (eq.(2.15)).
Thus, the variational formulation is equivalent to the weighted mean formulation with
linearized H. Since 3DVAR uses the time-independent background error covariance B ob-
tained from a long-term statistical mean, it is equivalent to OI as long as the observational
operator H is linear.
27
2.3 Data assimilation experiments on a simple system
In this section, we apply four data assimilation methods (3DVAR, serial EnSRF, LEKF,
full KF) on the Lorenz-96 model.
2.3.1 The Lorenz-96 model
The Lorenz-96 model (Lorenz 1996; Lorenz and Emanuel 1998) is deflned as
dxi
dt = xi?1(xi+1 ?xi?2)?xi +F (2.79)
Here, i = 1,???,N, where the boundary is cyclic, i.e. x?1 = xN?1, x0 = xN, and
x1 = xN+1. This model behaves chaotically in the case of external forcing F = 8.0,
in which case a time increment of 0.2 non-dimensional units corresponds to about one
day in terms of error growth rate (Lorenz 1996). The flrst term of the right hand side
simulates advection, and this model can be regarded as the time evolution of an arbitrary
one-dimensional quantity on a constant latitude circle, that is, the subscript i corresponds
to longitude. As in Lorenz (1996), we choose N = 40 and F = 8. The time step is chosen
as ?t = 0.01, and the time difierencing is solved using a 4th-order Runge-Kutta scheme.
2.3.2 Implementation of data assimilation methods
With a linear observational operator, 3DVAR and OI are equivalent, that is, they treat the
same problem with difierent approach as described in Section 2.2.8. Since the Lorenz-96
model has just 40 degrees of freedom, the full (40 ? 40) background covariance matrix
can be computed explicitly. Thus, in the simple model, 3DVAR and OI have the same
form, and the optimal weight (eq.(2.15)) is computed explicitly. To obtain the background
error covariance B, we apply the NMC method (Parrish and Derber 1992), accumulating
28
difierences between 24-hour forecasts and 18-hour forecasts valid at the same time as
samples of background errors.
In this simple system, the flve equations for Kalman flltering can be solved explic-
itly. Thus, the ?full Kalman flltering? (full-KF) without low-rank assumption has been
implemented.
As for the serial EnSRF and LEKF, model-independent core modules have been
developed. The details are described in Appendix A. Thus, only interfaces need to be
developed according to the model.
The implementation of serial EnSRF is straightforward. The system dimension in
core modules is set to 40. The core subroutines are called in the loop of observations
because observations are assimilated serially. After assimilating each observation, the
background ensemble is updated serially. After assimilating all observations, we retain the
flnal analysis ensemble. In serial EnSRF, we apply the multiplicative covariance in ation
described by eq.(2.66).
As for the implementation of LEKF, the local patch is taken as in Fig.2.2. We
compute the data assimilation equations for the larger region denoted by l. Since local
patches overlap each other, we average the overlapping values inside the smaller region
denoted by l2 to obtain the global analysis. In LEKF, we apply both the multiplicative
covariance in ation and the enhanced variance in ation (eq.(2.67)).
2.3.3 Description of experiments
The ?nature? or true run is created after long-term integration to avoid the initial spin-up
period. All of the data assimilation methods (3DVAR, full-KF, serial EnSRF, and LEKF)
are tested in the cases of 20 and 40 observations. We take the origin of the phase space
29
l2l
C
Figure 2.2: Deflnition of the local regions around a center point (C) for LEKF in the
Lorenz-96 model. The larger region denoted by l is the local patch of LEKF, the smaller
region denoted by l2 is the subset of the local patch that is averaged to obtain the global
analysis. Since l denotes the length of one side, the size of the local patch becomes 2l+1,
similarly for l2.
as the initial condition for the analysis experiments. The flrst 60 days in the analysis
experiments are considered as the spin-up period. To obtain the 1-year (365-day) mean
of the analysis root mean square error (RMSE), 425-day forecast-analysis cycle has been
performed. Observations are taken and analyzed every 6 hours. The observational error
standard deviation is 1.0 in the present experiments. The initial ensemble for serial EnSRF
and LEKF is chosen as vectors consisting of random numbers with Gaussian distribution.
Without data assimilation, the natural variability of the model shows an RMSE of 6.70,
which deflnes the error saturation limit.
2.3.4 Results in the case of 40 observations
Table 2.1 shows summary of the 1-year mean analysis RMSE in the case of 40 observations,
that is, all grid points have observations. All methods produce analysis state much closer
to the truth than observations. Observational error standard deviation is 1.0, about 15%
of the natural variability. All the Kalman flltering methods show the same performance
in their optimal cases (RMSE=0.2), whereas 3DVAR has errors with magnitudes about
30
3DVAR Full-KF Serial EnSRF LEKF
0.40 0.20 0.20 0.20
(? = 0.05) (? = 0.02,  = 5) (? = 0.015, l = 6, l2 = 2)
Table 2.1: The 1-year mean analysis RMSE of each method with optimal parameter
settings on the Lorenz-96 model when the number of observations is 40. Observational
error is given as 1.0, and the model natural variability without data assimilation shows
RMSE of 6.70. Ensemble size for the EnSRF and LEKF is flxed at 10.
twice as large (RMSE=0.40).
When serial EnSRF uses 40 ensemble members, the fllter is stable even without
the Schur product to localize the Kalman gain. When ensemble size is reduced, the
Schur product is required. There are two parameters for serial EnSRF: the localization
length scale  and the covariance in ation factor ?. The dependence on the parameters is
shown graphically in Fig.2.3, which is similar to Fig.3(b) of Whitaker and Hamill (2002).
The localization length scale  is difierent from the covariance fllter length scale L in
Whitaker and Hamill (2002). Whitaker and Hamill (2002) deflned the length scale L so
the covariance becomes 0 when d > L where d denotes distance. Our deflnition is given
by eqs.(2.49) and (2.50), where the distance d that yields zero covariance is given by
d > 2?
r10
3 ? (2.80)
Thus, the relationship between the length scale L in Whitaker and Hamill (2002) and the
length scale  used in this experiments is given as
L = 2?
r10
3 ? 
? 3.65 (2.81)
31
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Schur product length scale
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=40 Serial EnKF (NBV=10)
0.21
0.21
0.225
0.225
0.225
0.225
0.25
0.25
0.25
0.25
0.275
0.275
0.275
0.3
0.3
0.3
0.325
FILTER DIVERGENCE 
Figure 2.3: Parameter dependence of the analysis RMSE of serial EnSRF with 10 ensemble
members on the Lorenz-96 model when the number of observations is 40, cf. Fig.3(b) of
Whitaker and Hamill (2002). The horizontal and vertical axes show the localization length
scale  and the covariance in ation paramter ?, respectively. The minimum error 0.20 is
observed when  = 5 and ? = 0.02. ?FILTER DIVERGENCE? denotes the region with
RMSE of more than 1.0.
In considering the difierent measure of the length scale, Fig.2.3 shows almost identical
performance as Fig.3(b) in Whitaker and Hamill (2002), which supports the correctness
of the Fortran code developed here for the present research. According to Whitaker and
Hamill (2002), the minimum error was 0.20 which is identical to our results.
Fig.2.4 shows a similar flgure as Fig.2.3 in the case of LEKF. LEKF has three param-
eters: the local patch size 2l +1, the size of the local patch averaged to obtain the global
analysis 2l2 + 1, and the enhanced variance in ation factor ?. Since no clear dependence
32
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Size of local patch (l)
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=40 LEKF (NBV=10)
0.21
0.21
0.225
0.225
0.225
0.25
0.25
0.25
0.25
0.25
0.275
0.275
0.2750.3
FILTER DIVERGENCE 
Figure 2.4: The similar as Fig.2.3 but using LEKF with 10 ensemble members using
the enhanced variance in ation. The horizontal and vertical axes show the local patch
parameter l and the enhanced variance in ation paramter ?, respectively. l2 is flxed to 1.
The minimum error 0.20 is observed when l = 5 and ? = 0.01.
has been observed with respect to l2, l2 is flxed to 1. Note that the direct comparison
between serial EnSRF and LEKF may be misinterpreted since the way of localization is
difierent between serial EnSRF and LEKF, so that the metric in the horizontal axis may
be difierent. The minimum error 0.20 is observed when l = 5 and ? = 0.01.
In both EnKF, the following characteristics are observed:
? The smaller the localization length scale, the more stable the fllter (i.e., the farther
away from the region of fllter divergence).
? The larger the covariance in ation factor, the more stable the fllter.
33
3DVAR Full-KF Serial EnSRF LEKF
1.15 0.33 0.33 0.33
(B = 2.0?B) (? = 0.1) (? = 0.04,  = 6) (? = 0.025, l = 6, l2 = 2)
Table 2.2: The same as Table 2.1 but in the case of 20 observations and 10 ensemble
members.
? A larger length scale requires a larger covariance in ation factor for stable flltering.
? Filter divergence occurs with larger length scales and smaller in ation factors.
? The minimum error is observed in a medium length scale (? 5,6) with a smaller
in ation factor (? 0.02) near the boundary of the fllter divergence area.
2.3.5 Results in the case of 20 observations
The results in the case of 20 observations are summarized in Table 2.2 when the ensemble
size for EnKF is flxed to 10. The 20 observations are located regularly every other point.
Since the error variance is larger than in the case of 40 observations, the background
error covariance used in 3DVAR was multiplied by 2.0, which is a good choice among
several values we tried, though this may not be strictly optimal. The RMSE of 3DVAR is
larger than the observational error, whereas those of Kalman fllters are much smaller. The
difierence between 3DVAR and Kalman fllters becomes larger here, thus, Kalman fllters
show a larger advantage with fewer observations. All the Kalman flltering methods show
the same performances (RMSE=0.33) with their optimal settings.
Fig.2.5 shows the same flgure as Fig.2.3 but in the case of 20 observations. Compared
to the case of 40 observations, the fllter becomes considerably more unstable and the region
?FILTER DIVERGENCE? becomes larger, although the basic characteristic is similar.
34
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Schur product length scale
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=10)
0.35
0.35
0.375
0.375
0.375
0.375
0.4
0.4
0.4
0.425
0.425
0.425
0.45
0.45
0.45
0.475
0.475
0.475
0.5
0.5
0.5
FILTER DIVERGENCE 
Figure 2.5: The same as Fig.2.3 but in the case of 20 observations and serial EnSRF with
10 ensemble members. The minimum error 0.33 is observed when  = 6 and ? = 0.04.
The minimum error 0.33 occurs when  = 6 and ? = 0.04.
Fig.2.6 shows the same flgure as Fig.2.5 but ensemble size is increased to 20. When
ensemble size is increased, the fllter becomes more stable. In addition, the error becomes
slightly smaller almost everywhere. The minimum error 0.30 occurs when  = 6 and
? = 0.02.
Both serial EnSRF and LEKF work with 8 ensemble members. Fig.2.7 shows the
same flgure as Fig.2.5 but with 8 ensemble members. Fig.2.7 shows a similar structure as
Fig.2.5 but with larger ?fllter divergence? area and with larger errors almost everywhere.
Smaller ensemble size requires smaller covariance length scales. The minimum error 0.34 is
observed when  = 5 and ? = 0.04, which is slightly worse than in the case of 10 ensemble
members.
35
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Schur product length scale
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=20)
0.325
0.325
0.325
0.35
0.35
0.35
0.375
0.375
0.375 0.375
0.4
0.4
0.4
0.4
0.425
0.425
0.425
0.425
0.45
0.45
0.45
0.475
0.475
0.5
0.5
FILTER DIVERGENCE 
Figure 2.6: The same as Fig.2.5 but serial EnSRF with 20 ensemble members. The
minimum error 0.30 is observed when  = 6 and ? = 0.02.
36
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Schur product length scale
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=8)
0.35
0.350.375
0.375
0.4
0.4
0.4
0.40.425
0.425
0.4250.45
0.45
0.45
0.475
0.475
0.475
0.5
FILTER DIVERGENCE 
Figure 2.7: The same as Fig.2.5 but serial EnSRF with 8 ensemble members. The mini-
mum error 0.34 is observed when  = 5 and ? = 0.04.
37
l2 / in ation 0.03 0.04 0.05 0.06 0.07
1 .41 .37 .35 .36 .38
2 .34 .34 .35 .35 .37
3 .82 .34 .34 .38 .35
4 .34 .35 .40 .35 .36
5 .33 .35 .33 .35 .35
6 .81 .35 .35 .38 .36
Table 2.3: Analysis RMSE using LEKF with 8 ensemble members on Lorenz-96 model
when the number of observations is 20. Enhanced variance in ation is applied, and the
local patch parameter is flxed as l = 6. The other two parameters, l2 and ?, are changed.
Fig.2.8 shows the same flgure in the case of LEKF. Here, l2 is flxed to 0. It seems
LEKF has a smaller ?fllter divergence? area, thus, LEKF seems more stable than serial
EnSRF. However, the errors are slightly larger in LEKF. The minimum error 0.35 occurs
when l = 5 and ? = 0.03, which is slightly worse than serial EnSRF. As shown in Table
2.3, when we change l2 and ? with l flxed as l = 6, we get the minimum error 0.33 in
LEKF.
Although a smaller error appeared by changing l2, no clear dependence with respect
to l2 can be seen in Table 2.3. The  uctuation of the values suggests that the lengths of the
runs may not be su?cient to get reliable RMSE values, although the same computations
with difierent random numbers for observational errors performed several times yielded
the same overall results. Table 2.3 suggests that the LEKF errors are rather insensitive to
l2, since the only two large RMSE values (about 0.8) occur with the in ation parameter
? = 0.03 close to the area in which the fllter becomes divergent (Fig.2.8).
38
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Size of local patch (l)
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8)
0.375
0.4
0.40.4
0.4
0.425
0.425
0.425
0.425
0.425
0.45
0.45
0.45
0.475
0.475
0.475
0.5
0.5
0.5
0.525
FILTER DIVERGENCE 
Figure 2.8: The same as Fig.2.4 but LEKF with 8 ensemble members and enhanced
variance in ation. l2 is flxed to 0. The minimum error 0.35 is observed when l = 5 and
? = 0.03.
39
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Size of local patch (l)
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8)
0.4
0.425
0.425
0.425
0.425
0.4250.45
0.45
0.45
0.45
0.475
0.475
0.475
0.475
0.5
0.5
0.5
0.525
0.55
FILTER DIVERGENCE 
Figure 2.9: The same as Fig.2.8 but LEKF with 8 ensemble members and multiplicative
variance in ation. The way to in ate covariance is the same as that of serial EnSRF. The
minimum error 0.37 is observed when  = 5 and ? = 0.05.
Fig.2.9 shows the same flgure as Fig.2.8 but without enhanced variance in ation,
replaced by the same multiplicative covariance in ation used for the serial EnSRF. En-
hanced variance in ation contributes to stabilize the fllter, the minimum error becomes
smaller with the smaller in ation factor. In addition, with enhanced variance in ation, the
errors are smaller than with multiplicative in ation almost everywhere. Thus, enhanced
variance in ation contributes positively to the fllter performance.
We also applied observational error covariance localization within LEKF. Fig.2.10
shows the same flgure as Fig.2.9 but with observational error covariance localization.
Here, we flxed local patch parameters to l = 15 and l2 = 0, and varied localization scale
 of the observational error covariance localization. With observational error covariance
40
1 2 3 4 5 6 7 8 9 100.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Schur product length scale
covariance inflation factor
Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8) w/LOCALIZATION
0.35
0.35
0.35
0.375
0.375
0.375
0.375
0.4
0.4
0.4
0.4
0.425
0.425
0.425
0.45
0.45
0.450.475
0.475
0.5
0.5
FILTER DIVERGENCE 
Figure 2.10: The same as Fig.2.9 but LEKF with 8 ensemble members, multiplicative
variance in ation, and observational error covariance localization. The way to in ate
covariance is the same as that of serial EnSRF. The local patch parameters are flxed to
l = 15 and l2 = 0. The horizontal axis shows the length scale  of observational error
covariance localization. The minimum error 0.33 is observed when  = 4 and ? = 0.03.
localization, we obtain the lowest minimum error of 0.33 among all methods in the case of 8
ensemble members. Comparing Figs.2.10 and 2.7, we see that the use of localization allows
LEKF to use larger patches, with more observations, multiplicative in ation (instead of
enhanced in ation), and an overall performance at least as good as that of the serial
EnSRF with similar parameters. This observation has important implications for the
practical implementation of LEKF.
41
2.3.6 Online estimation of covariance in ation parameter
Kalnay (pers. comm.) suggested the covariance in ation parameter can be estimated by
the following statistics:
D
dd>
E
= (1+?)HPfH> +R (2.82)
a similar equation is found in eq.(3) of Houtekamer et al. (2005). Here, d denotes the
observational increment deflned as
d = yo ?Hxf (2.83)
Eq.(2.82) shows that the observational increment is the combination of the observational
error and the forecast error. To get the best estimate of the in ation parameter ?, one
can run the data assimilation cycle once and compute the statistics eq.(2.82) to update ?.
As an analogy of Kalman flltering in a model output statistics (MOS) (see for example,
appendix C of Kalnay 2003), we can apply Kalman flltering to make an online estimation
of ?. At each analysis time step, we estimate ? from eq.(2.82) to be ?o. We assume error
variance  o in this estimate, and we use the estimate at the previous time as a background
?b with error variance  b. Then, the analysis equations are given as:
?a = ?
b o +?o b
 b + o (2.84)
 a =
?
1?  
b
 b + o
!
 b (2.85)
This analyzed ?a can be used for the variance in ation parameter at the current step. ?a
and  a are used as background at the next time step.
We applied online estimation of covariance in ation parameter. Table 2.4 shows
analysis RMSE in the case of 20 observations using LEKF with 8 ensemble members.
Multiplicative covariance in ation and enhanced variance in ation are applied in both
42
local patch (l) or localization scale ( ) 1 2 3 4 5 6 7
Multiplicative .50 .46 .38 .38 .35 .36 .37
Enhanced .47 .40 .36 .37 .35 .35 .38
Multiplicative w/obs error localization .47 .38 .37 .34 .33 .34 .36
Enhanced w/obs error localization .49 .41 .36 .34 .34 .34 .38
Table 2.4: Analysis RMSE using LEKF with 8 ensemble members on Lorenz-96 model
when the number of observations is 20. Online estimation of variance in ation parameter is
applied. l2 is flxed to l2 = 0. In the cases with observational error covariance localization,
l is flxed to l = 15.
cases with and without observational error covariance localization. The online estimation
provides quite stable performance, although the o?ine optimally chosen constant in ation
parameter provides similar results. The smallest overall analysis error 0.33 is observed
with multiplicative in ation, a local patch parameter of l = 15, and an observational error
covariance localization size of 5.
This result is important because it provides a new approach to estimate covariance
in ation adaptively: In a global model the use of the same in ation for both tropical and
midlatitude regions is probably suboptimal since forecast errors (including model errors)
are larger in the tropics (e.g., Szunyogh et al, 2005).
2.3.7 Timing results
Computational time on a Linux PC with a 2.7GHz Intel Celeron (Northwood) processor is
shown in Table 2.5. The computation is a whole 425-day forecast-analysis cycle, containing
1700 data assimilations. The only difierence between 3DVAR and full-KF is covariance
43
3DVAR Full-KF Serial EnSRF LEKF
2.9 3.1 0.7 8.7
Table 2.5: Computational time in seconds on a Linux PC with an Intel Celeron 2.7GHz
processor for a 425-day forecast-analysis cycle on the Lorenz-96 model.
forecasts. Serial EnSRF is faster because it does not require matrix inversion. LEKF is
slower than full KF for the Lorenz-96 system.
Note that this timing result is very unrepresentative of realistic atmospheric systems
with high dimensionality, large number of observations and massively parallel computers.
LEKF has an advantage for parallel computation since each local patch is completely in-
dependent, the LEKF computation can be parallelized very e?ciently. In addition, LEKF
can treat large number of observations with less additional cost. Furthermore, the new lo-
cal ensemble transform KF (LETKF, Hunt 2005) avoids the eigenvalue decomposition and
is several times faster than LEKF (Harlim, Szunyogh, pers. comm.). On the other hand,
serial EnSRF cannot be parallelized e?ciently. Moreover, serial EnSRF may be expensive
in practice since it contains a loop over the number of observations. Implementation of
3DVAR is completely difierent for large systems, so the timing cannot be generalized.
2.3.8 Summary
The results on the Lorenz-96 model can be summarized as follows:
? The core modules of both serial EnSRF and LEKF are correctly coded since we have
been able to reproduce the results of Whitaker and Hamill (2002) and of Ott et al
(2004).
? All the Kalman flltering methods (full-KF, serial EnSRF, and LEKF) outperform
44
3DVAR, and the advantage becomes larger with fewer observations.
? All the Kalman flltering methods (full-KF, serial EnSRF, and LEKF) have compa-
rable performance with their optimal settings.
? Observational error covariance localization contributes to better LEKF performance
with little additional computational cost.
? Enhanced variance in ation in LEKF contributes to better fllter performance, but
multiplicative in ation becomes better when using localization.
? With the same multiplicative variance in ation, serial EnSRF outperforms LEKF,
but when localization and larger patches are used, LEKF is equal or better than
serial EnSRF with similar parameters.
? Smaller localization length scales and larger covariance in ation factors stabilize
EnKF, although both give suboptimal flltering.
? A new method introduced to estimate adaptively the optimal in ation factor works
well and gives results similar to the best obtained with o?ine tuning of a constant
in ation. In global models it could be used to estimate adaptive in ation factors
depending on location and time.
45
Chapter 3
Data assimilation on the SPEEDY primitive-equation model
3.1 Introduction
3.1.1 Overview
In Chapter 2, we described theoretical review and applied 3DVAR, serial EnSRF, LEKF,
and full KF on the Lorenz-96 model. All KF methods (serial EnSRF, LEKF, and full
KF) show similar results on the Lorenz-96 model, outperforming 3DVAR. In this chapter,
we describe the implementation and numerical experiments of three data assimilation
methods (3DVAR, serial EnSRF, and LEKF) on the SPEEDY primitive-equation model.
We explore the following questions in this chapter:
1. How does each method work on the SPEEDY model?
2. WhatistherelativeadvantagesanddisadvantagesbetweenserialEnSRFandLEKF?
3. How sensitive is the data assimilation to experimental settings?
4. What are the characteristics of the analysis and forecast error flelds?
A goal of this chapter is to investigate and compare the performance of the three
methods. In addition, we perform sensitivity experiments with difierent observational
networks, moisture observations, vertical error correlations, error covariance localization,
and random perturbation addition to ensemble members. These sensitivity experiments
investigate robustness of our results to the experimental settings. Moreover, they suggest
possible ways to improve the fllter performance. We also see the characteristics of the
46
analysis and forecast error flelds, suggesting how EnKF works.
3.1.2 The SPEEDY model
The SPEEDY model (Molteni 2003) is a recently developed atmospheric general circula-
tion model (AGCM) with a spectral primitive-equation dynamic core and a set of sim-
plifled physical parameterization schemes (SPEEDY stands for Simplifled Parameteriza-
tions, primitivE-Equation DYnamics). The goal of this model is to achieve computational
e?ciency while maintaining characteristics similar to the state-of-the-art AGCMs with
complex physics. The resolution of the model is T30L7 (horizontal spectral truncation of
30 wave numbers and 7 vertical levels), the computational cost is one order of magnitude
less than that of state-of-the-art AGCMs at similar horizontal resolution. According to
Molteni (2003), the SPEEDY model simulates the general structure of global atmospheric
circulation fairly well, and some aspects of the systematic errors are similar to many
AGCMs, though the error amplitude is larger than state-of-the-art models.
The SPEEDY model includes basic components of physical parameterizations used
in more complex GCMs, such as convection (a simplifled mass- ux scheme), large-scale
condensation, clouds, short-wave radiation (two spectral bands), long wave radiation (four
spectral bands), surface  uxes of momentum and energy (bulk aerodynamic formula),
and vertical difiusion. Details of the simplifled physical parameterization schemes of the
SPEEDY model can be found in Molteni (2003), especially in its Appendix which is
available on the website: ?http://www.ictp.trieste.it/?moltenif/speedy-doc.html?. The
boundary conditions of the SPEEDY model includes topographic height and land-sea
mask, which are constant, and sea surface temperature (SST), sea ice fraction, surface
temperature in the top soil layer, moisture in the top soil layer and the root-zone layer,
47
and snow depth, all of which are specifled by monthly means, and bare-surface albedo
and fraction of land-surface covered by vegetation, which are specifled by annual-mean
flelds. The lower boundary conditions such as SST are obtained by ECMWF?s reanalysis
in the period 1981-90. The incoming solar radiation  ux and the boundary conditions
(SST etc.), except bare-surface albedo and vegetation fraction, are updated daily.
Since the SPEEDY model is designed for long-term climate variability studies that
require ensembles of long-term integrations, the original code contains only time-mean
outputs. Thus, the input/output processes of a grid point value at intermittent time steps
had to be developed in order to enable short-term integrations of the forecast-analysis
cycle. The prognostic variables are zonal and meridional wind velocity components (u, v),
temperature (T), speciflc humidity (q), and surface pressure (ps). The grid point value
output is in the physical space with the grid size of 96?48?7. Inputs are only taken in
sigma levels, but outputs are in both sigma levels and pressure levels. The speciflc heights
of the vertical levels are shown in Table 3.1.
As for the timing results of running the model, a 3-month integration took about 6
minutes on a Linux PC with a 2.7GHz Intel Celeron (Northwood) processor, whereas it
took around 8 minutes for 24-hour integration of the Japanese operational global weather
prediction model with a reduced resolution (T42) on the same PC. For 6-hour cycle ex-
periments, a 6-hour forecast of the SPEEDY model requires about 2 seconds on the same
PC. Thus, a 2-month cycle experiment requires about 8-minute forecast computations for
each ensemble member.
48
Level index Sigma heights ( ) Pressure heights (hPa)
1 0.950 925
2 0.835 850
3 0.685 700
4 0.510 500
5 0.340 300
6 0.200 200
7 0.080 100
Table 3.1: Vertical levels of the SPEEDY model outputs. Sigma levels are also used as
model levels.
3.2 3DVAR implementation
3.2.1 Theory for practical implementation
The variational formulation has been introduced in Section 2.2.8. Practically, it is impos-
sible to implement the variational formulation explicitly because of the large number of
degrees of freedom of the system. Usually, NWP models have prognostic variables with
dimensions of O(107), thus, the matrix B has O(1014) elements, which requires at least
? 10TB of memory just for storing the matrix. However, most components of the matrix
B are very close to 0, and we can simplify B very much under reasonable assumptions.
We can understand this fact by considering the forecast error of temperature at College
Park, for example, has nothing to do with the forecast error of wind in Tokyo.
Parrish and Derber (1992) assumed zero spatial error correlation in the spectral
space. With the spectral transformation, they reduced most components of B. Alter-
49
natively, Barker et al. (2004) did not use spectral transformation in their 3DVAR built
for the nonhydrostatic flfth-generation Pennsylvania State University/National Center for
Atmospheric Research Mesoscale Model (MM5). We follow a similar approach as Barker
et al.
We deflne a variable transformation U as follows:
?x = U?v (3.1)
Here, ?xis deflned by eq.(2.75), ?v is a vector with the same dimension as ?x. Substituting
eq.(3.1) into eq.(2.73), we get
J(?v) = 12?v>U>B?1U?v+ 12(HU?v?d)>R?1(HU?v?d) (3.2)
In order for the covariance matrix to be the identity,
B = UU> (3.3)
should be satisfled. Then, eq.(3.2) can be rewritten as
J(?v) = 12?v>?v+ 12(HU?v?d)>R?1(HU?v?d) (3.4)
Difierencing eq.(3.4) with respect to ?v, we get the gradient of the cost function
rJ(?v) = ?v+U>H>R?1(HU?v?d) (3.5)
At this point, we can solve the 3DVAR problem for ?v. Since we start the minimizing
process from the background state, initially ?v = 0. Given the cost function (eq.(3.4)) and
its gradient (eq.(3.5)), a quasi-Newton minimizer flnds the solution for ?v, and eventually,
eq.(3.1) converts ?v to ?x, the analysis state in the real space. Note that only U and U>
are used in this process; U?1 is not required.
50
The construction of the variable transformation U is the key in the 3DVAR algo-
rithm. Usually it is assumed that U can be separated into a spatial correlation component
and an inter-variable correlation component, that is,
U = VCA (3.6)
where A, C and V stand for the error standard deviation, the spatial error correlation
and the inter-variable error correlation, respectively. Moreover, it is assumed C can be
separated into horizontal and vertical correlations. Usually, the horizontal correlation is
assumed to be a Gaussian shape. Only correlation length scale parameter needs to be
stored. Furthermore, the length scale is assumed to be constant in all points at the same
vertical level. Quasi-isotropy is also assumed. Full vertical correlation could be considered,
but it also can be assumed to be Gaussian.
An important part of 3DVAR is how to express the inter-variable correlation V,
which transforms prognostic variables into control variables that are assumed to be inde-
pendent of each other. Some prognostic variables are strongly dependent on each other
mostly because of the geostrophic balance. Usually in global models, only the geostrophic
balance is considered since it is the strongest and the most important balance in the
system.
3.2.2 Implementation on the SPEEDY model
In the previous section, we described a general idea to implement 3DVAR on a practical
atmospheric model following the approach by Barker et al. (2004). Here, we describe the
design of the variable transformation U in the present implementation.
The flrst part is the error standard deviationAin eq.(3.6). Since the present 3DVAR
analyzes not in the spectral space but in the physical grid space, we consider the full spatial
51
dependence of the error standard deviation. Almost no additional computation is required
by considering the full spatial dependence.
For the spatial error correlationCin eq.(3.6), the spatial correlation is separated into
vertical and horizontal correlations. The vertical resolution of the SPEEDY model is so
coarse that vertical background error correlation is not considered (experiments including
3 levels yielded worse results than one layer). Horizontal error correlation is assumed to
be Gaussian, and the 4th order recursive fllter (RF) technique (Purser et al. 2003) has
been implemented to estimate the Gaussian correlation shape. The 4th order RF is a
quasi-Gaussian fllter composed of a pair of forwarding and backwarding processes
Bi = flAi +fi1Bi?1 +fi2Bi?2 +fi3Bi?3 +fi4Bi?4 (3.7)
Ci = flBi +fi1Ci+1 +fi2Ci+2 +fi3Ci+3 +fi4Ci+4 (3.8)
where the subscript i denotes the space index of grid points in one direction. In eq.(3.7),
that is the forwarding process in the ascending order of i, A is input and the flltered value
B is obtained. In eq.(3.8), that is the backwarding process in the descending order of i, B
is input and the flltered value C is obtained. Eq.(3.7) and eq.(3.8) are mutually adjoint,
so eq.(3.7) may be used in U whereas eq.(3.8) may be used in U>. The coe?cients fl and
fi can be obtained by solving a 4th degree equation according to the correlation length
scale. The details of the recursive fllter technique are described in Appendix B.
For the inter-variable error correlation V in eq.(3.6), we consider the geostrophic
balance and introduce control variables as follows:
uu = u?r1ug(ps,T) (3.9)
vu = v?r2vg(ps,T) (3.10)
where r1 and r2 are regression coe?cients determined from statistics, and ug and vg are
52
zonal and meridional components of the geostrophic wind computed from ps and T using
the geostrophic balance equation on the sigma coordinate system (cf. Kalnay 2003, p.65,
eq.(2.6.27)):
fk?vg = ?RTrlnps ?r` (3.11)
where f, k, R, and ` denote Coriolis parameter, vertical unit vector, gas constant, and
geopotential height, respectively. Other variables (T, q, and ps) are assumed to be in-
dependent of the unbalanced wind (uu and vu). Thus, only the geostrophic balance is
considered for inter-variable correlation.
3.2.3 Background error statistics - NMC method
The design of variable transformation U described in the previous section requires back-
ground error statistics. We need to specify the error standard deviation, the length scale
of horizontal error correlations, and the regression coe?cients from the statistics. Back-
ground error statistics were computed using the NMC method (Parrish and Derber 1992),
accumulating the difierence between 24-hour forecast and 18-hour forecast as a sample of
the background error. In the present experiment, 87 samples from January 10 to January
31 were accumulated for the statistics. The analysis flelds for the flrst time are obtained
using 3DVAR with a reasonable guess of background error statistics, such as 500km for the
horizontal background error correlation length scale and the background error standard
deviation with the same magnitude as observational error standard deviation. The NMC
method provides better statistics. Once we get a more reasonable guess of background
error statistics using the NMC method, this process may be repeated.
Fig.3.1 shows the zonal mean of the regression coe?cients r1 with respect to u-
wind (eq.(3.9)). Almost no geostrophic balance is observed in the tropics, whereas strong
53
Figure 3.1: Zonal mean of the regression coe?cients r1 of u. The vertical and horizontal
axes show height and latitude respectively. The coe?cients re ect statistical strength of
the geostrophic relationship. The larger values, i.e. stronger geostrophic balance, are
observed in mid-latitudes.
geostrophy is observed in the mid-latitudes especially at the jet levels. The regression coef-
flcients r2 with respect to v-wind show a similar structure. Fig.3.2 shows the error standard
deviation of ps (top panel) and u at the 4th level (bottom panel), indicating strong spa-
tial dependence, especially in zonal structures in error standard deviation. Fig.3.3 shows
the vertical structure and latitudinal dependence of horizontal error correlations. There
is almost no vertical dependence (top panel), all levels except the top show almost no
correlation beyond 2 grid points. There is large latitudinal dependence in grid spacing
(bottom panel). This re ects the skewness of the global map, grid spacing is physically
more dense in higher latitudes.
The flgures shown as background error statistics so far are only small parts of the
54
Figure 3.2: Background error standard deviation of ps (top panel) and u at the 4th level
(bottom panel) that are estimated using the NMC method. The units for ps and u are Pa
and m/s, respectively. Strong spatial dependence especially in zonal structures is observed,
but their noisiness indicates sampling errors.
55
Figure 3.3: Background error horizontal correlations. The top panel shows vertical struc-
ture of u in y-direction, where the vertical axis shows sigma levels. The bottom panel
shows the latitudinal structure of ps in x-direction, where the vertical axis shows lati-
tudes. The horizontal axis shows horizontal length scale in grid spacing. The shades show
correlation.
56
background error information that are considered in the present 3DVAR. The present
3DVAR considers all components of standard deviation, the standard deviation has the
same size as a grid point value. As horizontal spatial correlation, zonal and meridional
components are stored separately. The zonal components are allowed to vary in latitude
and height, whereas meridional components are allowed to vary only in height. The
regression coe?cients measure the strength of the geostrophic balance and are allowed to
vary in height and latitudes. The regression coe?cients for u (r1) and v (r2) are separately
stored, although they should be equal in theory.
In order to examine our assumption of vertical independence, we compute vertical
correlations averaged horizontally. Fig.3.4 shows vertical correlations using the NMC
method. Some vertical correlations are seen at lower levels, but almost no correlation is
observed at upper levels. Except for moisture, only the flrst and second levels have large
correlations (about 0.4), other levels show less than 0.2 correlations. Thus, the vertical
independence is not a bad assumption in 3DVAR. Moisture, shown by long-short dashed
lines with black square, has largest vertical correlation at levels 2, 3, and 4, which may be
related to vertical moisture transport by large-scale convective parameterization.
3.2.4 Response tests
We perform a response test using a single observation of u with the observational increment
of 1.0m/s. Fig.3.5 shows responses when the observational station is located at mid-
latitude (top panel) and high-latitude (bottom panel) at the 4th level. It shows how
RF works, and it seems RF expands the observation signal to the Gaussian shape as
expected. In the higher latitude, the skewness of the map is considered correctly. Fig.3.6
shows analysis increments of wind vector and T (top panel) and ps (bottom panel) when
57
Figure 3.4: Vertical background error correlations of u-wind (solid lines with white circle),
v-wind(dashedlineswithblackcircle), temperature(shortdashedlineswithwhitesquare),
and moisture (long-short dashed lines with black square). The correlations are averaged
horizontally.
58
the observation is at mid-latitude at the 4th level. Although only u is observed, T and ps
are analyzed through the inter-variable correlation.
Since the single observation response tests show reasonable results, the 3DVAR is
tested using the dense observational network (Fig.3.9). Fig.3.7 shows analysis increment
(top panel) and the difierence between truth and flrst guess (bottom panel). Analysis
increment has a very similar structure as the difierence between truth and flrst guess,
which means the 3DVAR scheme succeeds as expected.
3.3 EnKF implementation
3.3.1 Serial EnSRF
We now apply the same core modules used for serial EnSRF in the Lorenz-96 model
(Section 2.3) to the SPEEDY model. Appendix A describes the core modules. The core
subroutines require one-dimensional arrays for inputs and outputs, all the variables are
combined together to form a one-dimensional array. All inter-variable correlations are
considered.
As described in Section 2.2.4, localization is very important to avoid sampling errors
caused by the limited ensemble size. As for horizontal localization, distance is computed
using the grid spacing unit. Thus, skewness of the map in the physical space is not consid-
ered. The boundary is cyclic in longitudinal direction, whereas in latitudinal direction, the
boundary is treated as a solid wall. The Schur product with the Gaussian-like weighting
function (eq.(2.49)) is applied in this grid spacing distance.
As in 3DVAR, errors in difierent vertical levels are considered to be independent.
Surface pressure (ps) is assumed to be correlated only with variables at the bottom level
( = 0.950). Thus, ps changes because of observations of ps itself and all other variables
59
Figure 3.5: Spatial response of the 3DVAR in the case of a single observation with the
observational increment of 1.0m/s at mid-latitude (top panel in m/s) and high-latitude
(bottom panel in m/s) at the 4th level ( = 0.51).
60
Figure 3.6: Inter-variable response of the 3DVAR in the case of a single observation with
the observational increment of 1.0m/s at mid-latitude at the 4th level. Analysis increment
of T (top panel in K) and ps (bottom panel in Pa) are shown.
61
Figure 3.7: Analysis increment of 3DVAR (top panel) and the difierence between truth
and flrst guess (bottom panel). The flelds are u-wind flelds (in m/s) at the 4th level.
62
l
l2
Figure 3.8: Deflnition of the local regions around a center point for LEKF in the SPEEDY
model. The larger region denoted by l is the local patch of LEKF, the smaller region
denoted by l2 is the subset of the local patch that is averaged to obtain the global analysis.
Since l denotes a half of a side, the area of the local patch becomes (2l + 1) ? (2l + 1),
similarly for l2.
(u, v, T, q) at the bottom level, and vice versa.
3.3.2 LEKF
In implementing the LEKF for the SPEEDY model, we use the same core modules from
the Lorenz-96 model (Section 2.3). Appendix A describes the core modules. The inputs
and outputs of the core subroutines are one-dimensional arrays. Thus, the model variables
are projected to the local space, and the local variables are combined together to be stored
in a one-dimensional array.
Fig.3.8 shows the local patch, which is a two-dimensional expansion of Fig.2.2 in
the case of the Lorenz-96 model. The region is a square in the grid space, the physical
shape is not a square because of the skewness of the map in higher latitudes. Since the
local patches are overlapped, only the values inside the smaller region denoted by l2 are
averaged to obtain the global analysis.
63
Cyclic boundary conditions are applied in the longitudinal direction. In the latitu-
dinal direction, the point next to a point (i,j) where i denotes longitude and j denotes
the highest latitude is taken as (i + 180,j). In this way, the same size of the local patch
is taken everywhere. This treatment is the same as Szunyogh et al. (2005).
As for the vertical localization, similar to the serial EnSRF, errors at difierent levels
are assumed to be independent except surface pressure (ps) unless otherwise noted. ps is
combined to other variables at all levels and thus, ps observations afiect other variables
at all levels. However, only the variables at the bottom level afiect ps. In this way, the
combination of ps and other variables is one-way in upper levels, whereas it is two-way
only at the bottom level. This treatment is also the same as Szunyogh et al. (2005).
In general, the vertical independence is not a good assumption, the system is de-
veloped so that it can consider vertical correlations. For comparison with 3DVAR and
serial EnSRF, vertical correlation is not considered, but the impact of vertical correlation
is shown later.
3.4 Experimental setup
So far, we described the implementation of three methods (3DVAR, serial EnSRF, and
LEKF). In this section, the condition of numerical experiments is described.
3.4.1 Observational network
We use three types of observational networks in the present research:
1. A regularly distributed dense network (Fig.3.9)
2. A regularly distributed sparse network (Fig.3.10)
3. An irregularly distributed realistic network (Fig.3.11)
64
Variable Observational error standard deviation
U-wind (u) 1.0 [m/s]
V-wind (v) 1.0 [m/s]
Temperature (T) 1.0 [K]
Speciflc humidity (q) 0.0001 [kg/kg]
Surface pressure (ps) 100.0 [Pa]
Table 3.2: Observational error standard deviation
Fig.3.9 shows the dense network that consists of 1056 total sounding stations (about 23
percent of all grid points) which are chosen regularly as 1 station out of every 4 grid
points excluding higher latitudes than 80 degrees north and south. Fig.3.10 shows the
sparse network that consists of 264 total sounding stations (about 6 percent of all grid
points) which are chosen regularly as 1 station out of every 16 grid points excluding higher
latitudes. Fig.3.11 shows the realistic observational network consists of 415 total sounding
stations (about 9 percent of all grid points) that simulate the real operational radiosonde
stations. The network is realistic in that more stations exist over land, especially in
Western Europe, Eastern Asia, and North America, with fewer stations over the ocean.
It is assumed that all sounding stations are located at model grid points so that no
interpolation is required. Each sounding station has complete vertical soundings, that is,
all variables are observed at all sigma levels. The observational error standard deviations
are shown in Table 3.2.
The dense observational network is dense enough for 3DVAR to work almost per-
fectly as shown later. Thus, in the dense network, there would be no signiflcant difierence
between EnKFs and 3DVAR. When the number of observations is reduced, advanced
65
Figure 3.9: A regularly distributed dense observational network. 1056 stations (23 percents
of all grid points) are located at one grid point out of every four grid points.
Figure 3.10: A regularly distributed sparse observational network. 264 stations (5.7 per-
cents of all grid points) are located at one grid point out of every 16 grid points.
66
Figure 3.11: An irregularly distributed realistic observational network. 415 stations (9.0
percents of all grid points) are located mostly over continents in the northern hemisphere.
methods are expected to show their advantage. In the sparse network, 3DVAR is ex-
pected to show large errors, whereas it is expected that EnKFs show much smaller errors.
Therefore, we compare the methods using the sparse network.
3.4.2 Data assimilation experiments
As described in Section 1.3, the data assimilation experiments are based on observational
systems simulation experiments (OSSEs) where the nature run is assumed to be known,
and observational data are simulated by adding random noise to the nature run according
to the observational errors. The nature run is created by the model integration after a
1-year spin-up period. By default, the SPEEDY model starts from the atmosphere at rest
(u = 0, v = 0) on the initial date January 1, 1981. We choose January 1, 1982 as the
initial date for data assimilation experiments after 1-year integration. We use the same
67
nature run for all data assimilation experiments in this chapter.
We generate another independent nature run starting from the atmosphere at rest
at a difierent time. The difierence between the two nature runs temporally averaged in
February, 1982 is 106.37 meters measured by the RMS difierence of 500hPa height fleld,
the error magnitude without data assimilation. We choose flelds at arbitrary times from
the second nature run to obtain the initial conditions of data assimilation. For 3DVAR,
a fleld is chosen. For EnKF, we choose flelds at arbitrary times for the initial ensemble.
The same initial ensemble is used for serial EnSRF and LEKF. If the ensemble size is
inflnitely large, the ensemble mean becomes the climatology. Thus, this is equivalent to
beginning the data assimilation cycle with climatology in 3DVAR. This choice is practical
in the sense that it is always available just by running the model, even if we have no idea
about the nature run.
Fig.3.12 shows a schematic of data assimilation experiments on the SPEEDY model.
The forecast fleld (?FORECAST?) by the SPEEDY model gives the flrst guess fleld for the
data assimilation system (?DA?), and the true fleld (?TRUTH?) with random noise pro-
vides observations. ?DA? may be 3DVAR, serial EnSRF, or LEKF. The output of ?DA?
is the analysis fleld (?ANALYSIS?) that provides the initial condition for the next model
integration. ?ANALYSIS? contains high frequency components, whereas ?INITIAL?, the
output by the SPEEDY model at forecast time 0, does not because the spectral transfor-
mation of the model smoothes out the high frequency components. Since the SPEEDY
model does not see the high frequency components, the smoothed ?INITIAL? fleld is con-
sidered as the analysis fleld for veriflcation. Thus, the analysis error is deflned as the
difierence between ?INITIAL? and ?TRUTH?.
We use the multiplicative covariance in ation (eq.(2.66)) both for serial EnSRF and
68
FORECAST
SPEEDY INITIAL
ANALYSIS
DA TRUTH
Add random noise as observational errors
Verification
Figure 3.12: Schematic of data assimilation experiments on the SPEEDY model. ?DA?
denotes the data assimilation system. The ?initial? fleld is the output by the SPEEDY
model which is difierent from the ?analysis? fleld because the high frequency components
are smoothed out by the spectral transformation of the model. The difierence between
?initial? and ?truth? is deflned as the analysis error.
69
LEKF. The enhanced variance in ation (eq.(2.67)) has been tested, but the fllter quickly
diverges in the SPEEDY model. The in ation factor is flxed to ? = 0.04, that is, 4%
covariance in ation, the same as in Szunyogh et al. (2005). This choice of in ation factor
is reasonable to make the EnKFs stable with reasonable choices of localization scale and
ensemble size.
3.5 Results
Following the development of the system, we perform analysis-forecast cycle experiments.
In this section we show the results of the three data assimilation systems (3DVAR, serial
EnSRF, LEKF) and compare them.
3.5.1 3DVAR
Fig.3.13 shows the analysis RMSE of 500hPa height fleld. The solid line shows the result
in the case of a dense observational network. The dense observational network is dense
enough for 3DVAR to perform almost perfectly; the analysis RMSE of 500hPa height is
less than 5 meters, much smaller than the operational NWP system for the real atmo-
sphere. RMSE averaged over one month after the initial one-month spin-up period is
3.58 meters. When the observations are decreased and the sparse network is applied, the
performance of 3DVAR worsens as expected. The dashed line shows the performance of
3DVAR using the same background error statistics optimal for the dense network. Since
the background error standard deviation of the sparse network is expected to be larger
than that of the dense network, the error standard deviation is multiplied by 2.0 giving
better performance (dotted line). The background error statistics have been regenerated
using the sparse network. The new background error information for the sparse network
70
Figure 3.13: Analysis RMSE of 500hPa height fleld in the case of 3DVAR with the dense
observational network (solid line), sparse observational network using the same background
error statistics based on the dense network (dashed line), sparse observational network
using background error standard deviation multiplied by 2.0 (dotted line), and sparse
observational network using background error statistics based on the sparse network (dot-
dashed line).
gives worse 3DVAR performance than the original background error information based
on the dense network (dot-dashed line). The three lines of 3DVAR with the background
statistics based on the dense network always outperform the run with no data assimilation,
which supports the usefulness of 3DVAR even in the sparse observational network.
The present 3DVAR system can consider full spatial dependencies of error standard
deviation, but we can easily simplify the dependencies by averaging the error standard
deviationinspace. Fig.3.14showsanalysisRMSEof500hPaheightinthecaseofthesparse
observational network with difierent spatial dependencies of the error standard deviation.
71
The solid line shows the results when we consider the full spatial dependencies of the error
standard deviation, the same as the dotted line in the previous Fig.3.13. The dashed
line shows the results in the case of horizontal mean of standard deviation, in which case
the standard deviation is a function of only the vertical levels and variables. The dotted
line shows the results in the case of zonal mean of standard deviation, in which case, the
standard deviation is a function of latitudes, vertical levels, and variables. All the three
performances are similar, although the horizontal mean gives better performance generally.
Table 3.3 shows the average performance of the three cases. The RMSE is averaged for one
month (112 samples) after the initial one-month spin-up period. The simpler structure of
background error standard deviation gives the better 3DVAR performance, suggesting that
the full spatial structure of the error covariance obtained with 87 samples is dominated
by sampling errors.
Since we observed the best performance (31.14 meters) when we multiplied the
background error standard deviation by 2.0 and took the horizontal mean of the original
error variance obtained from NMC-method statistics, hereafter, we use this conflguration
as an optimal setting of 3DVAR.
In order to see the error distributions in space, analysis RMSE of 500hPa height fleld
is plotted in space. Analysis errors are averaged only in time to obtain the spatial distri-
bution. Fig.3.15 shows the horizontal structure of analysis RMSE of 500hPa height fleld.
Shades and contours show analysis RMSE and mean analysis flelds, respectively. Dots
indicate observational locations. Large errors exist in polar regions where no observations
are available. We see a lattice-like pattern generated by the regularity of observational
locations. The top panel of Fig.3.16 shows the analysis error bias. We can also see large
analysis error bias in the polar region, and we can still flnd the lattice-like pattern. We
72
Figure 3.14: Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network with full spatial dependence of background error standard deviation (solid line),
with horizontally uniform background error standard deviation (dashed line), and with
zonally uniform background error standard deviation (dotted line).
3DVAR(full) 3DVAR(zonal mean) 3DVAR(hmean) NO ASSIM
RMSE[m] 36.65 34.62 31.14 106.37
Table 3.3: Analysis RMSE of 500hPa height fleld in meters using 3DVAR on the SPEEDY
model when the sparse observational network is applied. The RMSE is temporally aver-
aged for a month (112 samples) after the initial one-month spin-up period. ?3DVAR(full)?
denotes the case that full spatial dependence of background error standard deviation is
considered. ?3DVAR(zonal mean)? and ?3DVAR(hmean)? denote the cases that zonally
and horizontally uniform background error standard deviation is considered, respectively.
?NO ASSIM? denotes the case that no data assimilation is applied.
73
Figure 3.15: Horizontal structure of analysis RMSE of height fleld in the case of 3DVAR
with horizontally uniform background error standard deviation and the sparse observa-
tional network. Shades and contours show analysis RMSE (in meters) and mean analysis
flelds (in meters), respectively. Dots indicate observational locations.
cannot neglect the error bias in this one-month period even if the model is perfect. The
bias fleld shows a clear signal in a certain wavelengths, mostly the wavelength of distance
between observations. In the longitudinal direction in higher latitudes, longer wavelengths
show a stronger signal. It seems the regular observational network introduces a regular
bias structure. If we subtract error bias from the analysis error fleld and recalculate the
analysis RMSE, i.e. standard deviation, the lattice-like pattern of the analysis RMSE is
reduced as shown in the bottom panel of Fig.3.16. In addition, the large errors in polar
regions are reduced, thus, the error bias explains a large amount of the polar errors.
So far, we have examined the analysis RMSE of 500hPa height fleld. Although
this fleld is often used as a proxy of the performance of a forecast model and a data
74
Figure 3.16: Horizontal structure of analysis error bias (top panel) and standard deviation
(bottom panel) of height fleld (in meters) in the same case as Fig.3.15. Shades show
analysis error bias (top panel) and standard deviation (bottom panel), contours show
mean analysis flelds. Dots indicate observational locations.
75
assimilation system, it is important to see other variables at other levels. Fig.3.17 shows
the zonal structures of analysis RMSE of the u-wind (top panel) and temperature (bottom
panel). Shades and contours show analysis RMSE and mean analysis flelds, respectively.
If we ignore the large errors in the polar regions, we see large analysis error in the tropical
regions around the 200hPa level in both flgures. In addition, both flgures show striped
pattern, where the minima correpond to observational locations. For the temperature
fleld (bottom panel), we see the tendency for larger errors to appear in upper levels in the
tropics. The error is also large in lower levels in middle and high latitudes. This feature
is similar to what Szunyogh et al. (2005) obtained.
The results using 3DVAR are summarized as follows:
? The 3DVAR data assimilation cycle experiment using a dense observational network
shows almost perfect performance: 3.58m RMSE of 500hPa height fleld.
? The 3DVAR data assimilation cycle experiments using the sparse observational net-
work show much worse performance: 31.14m RMSE of 500hPa height fleld.
? Horizontally uniform background error standard deviation results in better 3DVAR
performance, although the 3DVAR system can consider full spatial dependence,
suggesting that a sample of 87 is not enough for the full spatial dependence.
? Background error statistics estimated with the dense observational network have a
good structure for data assimilation with the sparse observational network.
? Horizontal distribution of analysis RMSE shows a lattice-like pattern generated by
the regularity of observational locations.
? Zonal mean structure of analysis RMSE shows a striped pattern corresponding to
the observational locations.
76
Figure 3.17: Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of 3DVAR with horizontally uniform
background error standard deviation and the sparse observational network. Shades and
contours show analysis RMSE and mean analysis flelds, respectively.
77
NBV / Schur 1 2 3 4 5 6
10 11.74 24.60 40.17 D D D
20 10.80 10.69 10.99 D 22.37 D
30 15.07 5.47 4.99 8.00 10.48 17.02
Table 3.4: Analysis RMSE of 500hPa height fleld in meters using serial EnSRF on the
SPEEDY model when the sparse observational network is applied. The RMSE is tempo-
rally averaged for a month (112 samples) after the initial one-month spin-up period. For
comparison, the analysis RMSE of 3DVAR in the same period is 31.14. The ensemble size
(NBV) is chosen as 10, 20, or 30, and the horizontal length scale of the Schur product is
chosen as integers from 1 to 6. ?D? denotes fllter divergence. 4% multiplicative covariance
in ation is applied.
3.5.2 Serial EnSRF
Table 3.4 shows the analysis RMSE of 500hPa height fleld. The ensemble size is chosen
to be 10, 20, and 30, and the horizontal correlation length scale for the Schur product is
chosen from the integers between 1 and 6. Even with 10 ensemble members, serial EnSRF
outperforms 3DVAR whose best result shown in the previous section was 31.14m. The
best performance of serial EnSRF is seen in the case of 30 ensemble members with the
horizontal scale of 3.0, which shows almost perfect performance (as good as 3DVAR with
the dense observational network).
A smaller ensemble size requires a smaller horizontal length scale of the Schur prod-
uct. The data assimilation cycle is more stable for a larger ensemble size. Both results
are similar to those of the Lorenz-96 model.
Fig.3.18 shows temporal sequences of the analysis RMSE of 500hPa height fleld in
78
Figure 3.18: Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using serial EnSRF with 10 (solid line), 20 (dashed line), and 30 ensemble mem-
bers (dotted line). The localization length scale of the Schur product is chosen as optimal.
For comparison, the RMSE of 3DVAR is shown as the short dashed line.
the cases of serial EnSRF with 10 (solid line), 20 (dashed line), and 30 ensemble members
(dotted line). The localization length scale of the Schur product is chosen as optimal, that
gives the best result in each case in the Table 3.4 (marked by rectangles). For comparison,
the RMSE of 3DVAR is shown as the short dashed line. Serial EnSRF clearly outperfoms
3DVAR. The case of 30 ensemble members clearly outperforms others almost everywhere.
In order to see the error distributions in space, analysis RMSE of 500hPa height
fleld is plotted in space, corresponding to Figs.3.15 and 3.16 to be compared with 3DVAR.
Fig.3.19 shows the horizontal structure of analysis RMSE of 500hPa height fleld, similar
to Fig.3.15. Large errors exist in polar regions where no observation is available, but
the errors are smaller than those in 3DVAR. We can still see a lattice-like pattern cor-
79
Figure 3.19: Horizontal structure of analysis RMSE of height fleld, similar to Fig.3.15, in
the case of serial EnSRF with 20 ensemble members and localization scale of 2.0. Shades
show analysis RMSE (in meters), contours show mean analysis flelds (in meters). Dots
indicate observational locations.
responding to the observational locations. The top panel of Fig.3.20 shows the analysis
error bias. Basic characteristics remain the same as the case of 3DVAR: large errors with
low wavenumbers exist in polar regions. We flnd a lattice-like pattern with the same
wavelength as the distance between observational points. The bottom panel of Fig.3.20
shows analysis error standard deviation, it is obtained in the same manner as RMSE but
subtracting error bias. The standard deviation fleld shows almost no lattice-like pattern.
Thus, we can attribute the lattice-like pattern to the error bias. In addition, the large
errors in polar regions are reduced, thus, the error bias explains the majority of the polar
errors.
Fig.3.21 shows the analysis RMSE of four variables at all seven pressure levels. As in
80
Figure 3.20: Horizontal structure of analysis error bias (top panel) and standard deviation
(bottom panel) of height fleld (in meters), similar to Fig.3.16, in the case of serial EnSRF
with 20 ensemble members and localization scale of 2.0. Shades show analysis error bias
(top panel) and standard deviation (bottom panel), contours show mean analysis flelds.
Dots indicate observational locations.
81
Table 3.4, RMSE is averaged for one month after the initial one-month spin-up period (112
samples), so that we get a single number for one variable at one level. Fig.3.21(a) shows
the analysis RMSE of u-wind flelds. A smaller RMSE is observed at levels higher than
300hPa, this trend is clearer in 3DVAR. As expected, the larger ensemble size gives better
performance everywhere, although the cases of 10 and 20 ensemble members are similar.
The case of 30 ensemble members clearly outperforms the others. It is noted that only
the case of 30 ensemble members gives smaller RMSE than observational error standard
deviation (1.0m/s, shown by the thin solid line) at all levels. Fig.3.21(b) shows the analysis
RMSE of height flelds, the serial EnSRF clearly outperforms 3DVAR at all levels. Serial
EnSRF performs almost constantly at all levels, whereas 3DVAR performs better at upper
levels. Fig.3.21(c) shows the analysis RMSE of temperature flelds. Generally, higher
altitude has smaller RMSE, especially clear in 3DVAR. All cases of serial EnSRF provide
smaller RMSE than observational error standard deviation (1.0K, shown by the thin solid
line) at all levels. Fig.3.21(d) shows the analysis RMSE of speciflc humidity flelds. The
order of the performances is basically the same as the previous flgures (a)-(c). Higher
altitudes are drier, so the RMSE is smaller in higher altitudes. Still, data assimilation
provides useful information, and we can see the difierence among difierent data assimilation
schemes and settings in higher altitudes. Lower altitudes contain more humidity, and
mostly the RMSE exceeds the observational error standard deviation (10?4kg/kg, shown
by the thin solid line).
Fig.3.22 shows the zonal structures of analysis RMSE of the u-wind (top panel) and
temperature (bottom panel), corresponding to Fig.3.17 of 3DVAR. Shades and contours
show analysis RMSE and mean analysis flelds, respectively. It is noted that the errors in
the polar regions are much smaller than those in 3DVAR. Basic features seen in 3DVAR
82
Figure 3.21: Analysis RMSE at the all pressure levels temporally averaged for one month
(112 samples) after the initial one-month spin-up period in the case of the sparse observa-
tional network using serial EnSRF with 10 (solid line), 20 (dashed line), and 30 ensemble
members (dotted line). The four panels (a), (b), (c) and (d) correspond to u-wind, height,
temperature (T) and speciflc humidity (q) flelds, respectively. The observational error
standard deviations are shown as thin solid lines wherever applicable. The localization
length scale of the Schur product is chosen as optimal. For comparison, the RMSE of
3DVAR is shown as the short dashed line.
83
remain similar, but with smaller amplitudes. The striped pattern seen in 3DVAR also still
exists, but with much smaller amplitudes.
The results using serial EnSRF are summarized as follows:
? In the case of the sparse observational network, serial EnSRF clearly outperforms
3DVAR and shows almost perfect performance (4.99m RMSE of 500hPa height fleld
with an optimal setting versus 31.14 for 3DVAR).
? The larger ensemble size, the better the performance.
? A smaller ensemble size requires smaller localization length scale, as observed in the
Lorenz-96 model.
? Even with only 10 ensemble members, serial EnSRF outperforms 3DVAR.
? Horizontal and zonal distributions of analysis RMSE show lattice-like and striped
patterns similar to 3DVAR but with smaller amplitudes, we can attribute the lattice-
like pattern to the analysis error bias.
? Not only 500hPa height flelds, but also other flelds at other levels show much better
performance than 3DVAR.
3.5.3 LEKF
Table 3.5 shows the analysis RMSE of 500hPa height fleld, similar to Table 3.4 in the case
of serial EnSRF. The ensemble size is chosen to be 10, 20, and 30, and the local patch
parameter l is chosen from 1 to 6. LEKF outperforms 3DVAR (31.14m at the best) for
20 and 30 ensemble members. The best performance is seen in the case of 30 ensemble
members with the local patch parameter l = 4.
84
Figure 3.22: Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of serial EnSRF with 20 ensemble
members and localization scale of 2.0. Shades and contours show analysis RMSE and
mean analysis flelds, respectively.
85
NBV / l 1 2 3 4 5 6
10 60.38 75.88 D D D D
20 35.45 19.81 16.85 30.84 D D
30 17.68 17.40 8.67 8.44 13.23 14.50
Table 3.5: Analysis RMSE of 500hPa height fleld in meters using LEKF on the SPEEDY
model when the sparse observational network is applied. The RMSE is temporally aver-
aged for a month (112 samples) after the initial one-month spin-up period. For comparison,
the analysis RMSE of 3DVAR in the same period is 31.14. The ensemble size (NBV) is
chosen as 10, 20, or 30, and the local patch parameter l is chosen as integers from 1 to 6.
?D? denotes fllter divergence. 4% multiplicative covariance in ation is applied.
A smaller ensemble size requires a smaller local patch size. The data assimilation
cycle is more stable for larger ensemble size. Both results are the same as in the case of
the Lorenz-96 model.
Fig.3.23 shows an LEKF version of Fig.3.18, which shows temporal sequences of the
analysis RMSE of 500hPa height fleld in the cases of LEKF with 10 (solid line), 20 (dashed
line), and 30 ensemble members (dotted line). The local patch parameter l is chosen as
optimal, it gives the best result in each case in the Table 3.5 (surrounded by rectangles).
For comparison, the RMSE of 3DVAR is shown as the short dashed line. LEKF with 20
and 30 ensemble members clearly outperfoms 3DVAR. The case of 30 ensemble members
outperforms others almost everywhere.
In order to see the error distributions in space, analysis RMSE of 500hPa height
fleld is plotted in space, similar to Figs.3.15, 3.16, 3.19 and 3.20. Fig.3.24 shows the
horizontal structure of analysis RMSE of 500hPa height fleld, similar to Figs.3.15 and
86
Figure 3.23: Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using LEKF with 10 (solid line), 20 (dashed line), and 30 ensemble members
(dotted line). The local patch parameter l is chosen as optimal. For comparison, the
RMSE of 3DVAR is shown by the short dashed line.
87
3.19. Large errors exist in polar regions where no observations are available, but the
errors are slightly larger than those in serial EnSRF. We can see a lattice-like pattern
corresponding to the observational locations. The top panel of Fig.3.25 shows the analysis
error bias. Basic characteristics remain the same as the case of 3DVAR: large errors
with low wavenumbers exists in polar regions and we see a lattice-like pattern with the
same wavelength as distance between observational points. The bottom panel of Fig.3.25
shows analysis error standard deviation, which is obtained in the same way as obtaining
RMSE but subtracting error bias. The standard deviation fleld shows almost no lattice-
like pattern, as in the case of serial EnSRF. Thus, we can attribute the lattice-like pattern
to the error bias. In addition, the large errors in polar regions are reduced, thus, the error
bias explains large amount of the polar errors.
As was shown in the case of serial EnSRF (Fig.3.21), the RMSEs of other variables
at all levels in the case of LEKF are also shown here. LEKF with 10 ensemble members
is worse than 3DVAR everywhere, which is omitted from the following flgures. Thus,
3DVAR, LEKF with 20 and 30 ensemble members are shown by short dashed lines, solid
lines, and dashed lines, respectively.
Fig.3.26 shows the analysis RMSE of four variables at the all seven pressure levels.
Fig.3.26(a) shows the analysis RMSE of u-wind flelds. Smaller RMSE is observed at the
levels higher than 300hPa. This trend is clearer for 3DVAR. As expected, the larger
ensemble size performs better everywhere. It is noted that only the case of 30 ensemble
members gives smaller RMSE than observational error standard deviation (1.0m/s, shown
by the thin solid line) at the lowest three levels. LEKF with 30 ensemble members clearly
outperforms the others everywhere. Fig.3.26(b) shows the analysis RMSE of height flelds,
where we see that LEKF outperforms 3DVAR, especially at lower levels. At the top level,
88
Figure 3.24: Horizontal structure of analysis RMSE of height fleld, similar to Figs.3.15
and 3.19, in the case of LEKF with 20 ensemble members and local patch parameter l = 3.
Shades show analysis RMSE (in meters), contours show mean analysis flelds (in meters).
Dots indicate observational locations.
89
Figure 3.25: Horizontal structure of analysis error bias (top panel) and standard deviation
(bottom panel) of height fleld (in meters), similar to Figs.3.16 and 3.20, in the case of
LEKF with 20 ensemble members and local patch parameter l = 3. Shades show analysis
error bias (top panel) and standard deviation (bottom panel), contours show mean analysis
flelds. Dots indicate observational locations.
90
there is almost no difierence between 3DVAR and LEKF with 20 ensemble members.
LEKF with 30 ensemble members clearly outperforms the others everywhere. Fig.3.26(c)
shows the analysis RMSE of temperature flelds. Generally, higher altitude gives smaller
RMSE, that is clear in 3DVAR. Almost all cases of LEKF provide smaller RMSE than
observational error standard deviation (1.0K, shown by the thin solid line). LEKF with 30
ensemble members clearly outperforms the others everywhere, and shows almost constant
RMSE in the lower 6 levels. Fig.3.26(d) shows analysis RMSE of speciflc humidity flelds,
where we can see that higher altitude is drier, so the RMSE is smaller in higher altitudes.
All show similar performances, and no clear difierence can be seen, whereas there was clear
difierence in the case of serial EnSRF (cf. Fig.3.21(d)). Lower altitudes contains more
humidity, and the RMSEs exceed the observational error standard deviation (10?4kg/kg,
shown in the thin solid line).
Fig.3.27 shows the zonal structures of analysis RMSE of the u-wind (top panel) and
temperature (bottom panel), similar to Figs.3.17 and 3.22 of 3DVAR and serial EnSRF.
Shades and contours show analysis RMSE and mean analysis flelds, respectively. We can
see the same basic features seen in 3DVAR and serial EnSRF: large errors in u-wind fleld
exist in upper tropics, large errors in the temperature fleld exist in upper tropics and lower
mid-latitudes. The amplitude of the errors is smaller than 3DVAR, but larger than serial
EnSRF. The striped pattern seen in 3DVAR also still exists with much smaller amplitudes.
Table 3.6 shows the algorithm sensitivity to the parameter l2. The table is created
in the same way as Table 3.5, that is, the temporal mean of the analysis RMSE of 500hPa
height fleld is computed using various parameter settings. Here, the local patch parameter
is flxed to l = 2, and the ensemble size is flxed to 20 or 30. As in the Lorenz-96 model, we
see no clear dependence on the parameter l2, though the best performance in a given en-
91
Figure 3.26: Analysis RMSE at all pressure levels temporally averaged for one month (112
samples) after the initial one-month spin-up period in the case of the sparse observational
network using LEKF with 20 (solid line) and 30 ensemble members (dashed line). The
four panels (a), (b), (c) and (d) correspond to u-wind, height, temperature (T) and speciflc
humidity (q) flelds, respectively. The observational error standard deviations are shown
as thin solid lines wherever applicable. The localization length scale of the Schur product
is chosen as optimal. For comparison, the RMSE of 3DVAR is shown as the short dashed
line.
92
Figure 3.27: Zonal structure of analysis RMSE of u-wind fleld (top panel, in m/s) and
temperature fleld (bottom panel, in K) in the case of LEKF with 20 ensemble members
and local patch parameter l = 3. Shades and contours show analysis RMSE and mean
analysis flelds, respectively.
93
NBV / l2 0 1 2
20 19.81 19.22 17.24
30 17.40 23.81 10.70
Table 3.6: Analysis RMSE of 500hPa height fleld in meters using LEKF on the SPEEDY
model when the sparse observational network is applied and the parameter l2 is changed.
The ensemble size (NBV) is chosen to be 20 or 30, the local patch parameter is flxed
to l = 2. The RMSE is temporally averaged for a month (112 samples) after the initial
one-month spin-up period. For comparison, the analysis RMSE of 3DVAR in the same
period is 31.14.
semble size is observed at the largest l2 value (shown by rectangles). The large  uctuations
may be caused by the limited run length, although overall results look reasonable.
The results using LEKF are summarized as follows:
? In the case of the sparse observational network, LEKF clearly outperforms 3DVAR
and shows very good performance (8.44m RMSE of 500hPa height fleld with an
optimal setting).
? The larger ensemble size, the better the performance.
? A smaller ensemble size requires smaller localization length scale, as in the Lorenz-96
model.
? 10 ensemble members are not enough for the LEKF to outperform 3DVAR. How-
ever, with the observational error covariance localization, LEKF with 10 ensemble
members outperform 3DVAR.
? We observe no clear dependence on the parameter l2.
94
? Horizontal and zonal distributions of analysis RMSE show lattice-like and striped
patterns with smaller amplitudes than 3DVAR but with larger amplitudes than serial
EnSRF, we can attribute the lattice-like pattern to the analysis error bias.
? Not only 500hPa height flelds, but also most other flelds at other levels show better
performance than 3DVAR in the cases of 20 and 30 ensemble members.
3.5.4 Timing results
In 3DVAR, the computational time depends on the number of iterations of the variational
process. Each iteration requires less than a half second on a Linux PC with a 2.7GHz
Intel Celeron (Northwood) processor. A typical number of iterations until convergence is
about 35, which requires about 5 seconds of computational time.
Table 3.7 shows timing of serial EnSRF for each analysis. Here, the sparse observa-
tional network is used. The analysis is performed on two Linux PCs: one with a 2.7GHz
Intel Celeron (Northwood) and the other with a 2.53GHz Intel Celeron D (Prescott). The
latter is a later model with a twice larger L2-cache memory. Although the clock speed of
the former processor (2.7GHz) is larger, the latter one (2.53GHz) performs almost twice
as fast.
Computational time of serial EnSRF strongly depends on the number of observa-
tions and the ensemble size. However, it does not depend on the horizontal localization
scale. The system can be further accelerated by optimizing the code by considering the
localization scale.
The computational time of LEKF strongly depends on the ensemble size and the
local patch parameters(both l and l2). However, it does not strongly depend on the
number of observations. The system has been accelerated by optimizing the code as much
95
NBV 10 20 30
Celeron 2.7GHz (Northwood) 2:20 4:30
Celeron D 2.53GHz (Prescott) 1:15 2:20 3:20
Table 3.7: Timing (min:sec) of serial EnSRF with 10, 20, and 30 ensemble members on
the SPEEDY model when the sparse observational network is applied.
NBV(l) 10(1) 20(3) 30(4)
Celeron 2.7GHz (Northwood) 2:50 8:30
Celeron D 2.53GHz (Prescott) 1:35 6:00 16:30
Table 3.8: Timing (min:sec) of LEKF with 10, 20, and 30 ensemble members on the
SPEEDY model when the sparse observational network is applied. The local patch pa-
rameter l is chosen as an optimal setting shown in the parenthesis, and l2 is flxed to
0.
as possible. However, we do not apply the new technique to avoid eigenvalue decomposition
(local ensemble transform Kalman fllter, LETKF) proposed by Hunt (2005). Note that
LEKF can be parallelized very e?ciently, thus, LEKF can be accelerated almost perfectly
on parallel computers.
Table 3.8 shows timing of LEKF for each analysis. Here, the sparse observational
network is used, l2 is flxed to 0. The local patch parameter l is chosen as an optimal setting
that gives the best performance in a given ensemble size. The analysis is performed on
two Linux PCs: one with a 2.7GHz Celeron (Northwood) and the other with a 2.53GHz
Celeron D (Prescott), as in the case of serial EnSRF (Table 3.7). Similar to the case of
serial EnSRF, the latter processor performs faster.
96
Ensemble size (NBV) 10 20 30
Optimal localization scale  of serial EnSRF 1 2 3
Optimal local patch parameter l of LEKF 1 3 4
Table 3.9: Optimal localization parameters (length scale  of serial EnSRF and patch size
l of LEKF) with given ensemble sizes. Here, the covariance in ation parameter is flxed to
4%, the parameter l2 in LEKF is flxed to 0.
3.5.5 Comparison
For comparison of each method, the results of each method are combined and shown
in Table 3.10 and Fig.3.28. Each result has already been shown separately in previous
sections, the optimal parameters are chosen as shown in Table 3.9. Serial EnSRF is the
best among the three methods. This result is the same as that in the case of the Lorenz-96
model when the same multiplicative covariance in ation is applied in both serial EnSRF
and LEKF. Note that in the Lorenz-96 model, enhanced variance in ation signiflcantly
improved LEKF?s performance, LEKF is comparable to serial EnSRF. The RMSE values
contain errors due to the limited run length, but we obtained a similar value with doubled
run length in the case of serial EnSRF with 10 ensemble members. Fig.3.39 shows that
the rms errors of LEKF and EnSRF with 30 ensemble members are essentially identical
except poleward of 80 degrees, where the LEKF has larger errors which could be improved
with a difierent formulation for the polar patches in the LEKF.
As for computational time, Table 3.11 shows timing on a Linux PC with a 2.53GHz
Intel Celeron D (Prescott) processor. The timing was adopted from Tables 3.7 and 3.8. On
a single computer with serial treatment, it is clear that serial EnSRF is faster, especially
in the larger ensemble size.
97
NBV Serial EnSRF LEKF 3DVAR NO ASSIM
10 11.74 60.38 31.14 106.37
20 10.69 16.85 31.14 106.37
30 4.99 8.44 31.14 106.37
Table 3.10: Comparison of analysis RMSE of 500hPa height fleld in meters on the SPEEDY
model when the sparse observational network is applied. For each ensemble size (NBV),
best results are chosen for serial EnSRF and LEKF from Tables 3.4 and 3.5. For 3DVAR,
the best performance 31.14m is chosen. ?NO ASSIM? denotes the case without data
assimilation.
Figure 3.28: Analysis RMSE of 500hPa height fleld in the case of the sparse observational
network using 3DVAR (short dashed line), serial EnSRF (solid line) and LEKF (dashed
line). The ensemble size for serial EnSRF and LEKF is 30. Parameters are optimal, that
is, best results are shown for each method.
98
NBV 10 20 30
Serial EnSRF 1:15 2:20 3:20
LEKF 1:35 6:00 16:30
Table 3.11: Timing (min:sec) of serial EnSRF and LEKF with 10, 20, and 30 ensemble
members on the SPEEDY model when the sparse observational network is applied.
However, LEKF can be parallelized very e?ciently, thus, in parallel computers with
tens of computational nodes, LEKF could be much faster than serial EnSRF. Furthermore,
timing of serial EnSRF strongly depends on the number of observations unlike LEKF. More
satellite data will be available in the near future. With larger number of observations,
LEKF has an important advantage in terms of computational e?ciency.
Simulating the realistic situation, we estimate timing when 100 times more obser-
vations (26400 observations) are assimilated using a parallel computer with 100 computa-
tional nodes. We assume each node has the same performance as a 2.53GHz Intel Celeron
D (Prescott) processor. Let the number of computational nodes and parallelization e?-
ciency denote N and fi, respectively. The computation is accelerated as
T0 = T ? 1Nfi (3.12)
where T and T0 denote timings using 1 node and N nodes, respectively. LEKF is e?ciently
parallelized, we assume 80% e?ciency. Serial EnSRF could be parallelized by treating
distant observations separately by virtue of Schur product localization, we assume 20%
e?ciency. LEKF can treat large number of observations more e?ciently, we assume 50
times more computation caused by 100 times more observations. By contrast, Serial
EnSRF simply increases the number of loops, we assume 100 times more computation
caused by 100 times more observations. In summary, we multiply the timing results of
99
NBV 10 20 30
Serial EnSRF 6:15 11:40 16:40
LEKF 1:00 3:45 10:19
Table 3.12: Estimated timing (min:sec) of serial EnSRF and LEKF with 10, 20, and
30 ensemble members when 100 times more observations are assimilated using a parallel
computer with 100 computational nodes. Parallelization e?ciency is assumed as 80% and
20% for LEKF and serial EnSRF, respectively. We assume the increased observations
require 50 and 100 times more computation for LEKF and serial EnSRF, respectively.
Table 3.11 as follows:
T0LEKF = TLEKF ? 58 (3.13)
T0EnSRF = TEnSRF ?5 (3.14)
Table 3.12 shows the estimated timimg. By the e?cient parallelization and large number
of observations, LEKF is expected to be faster everywhere. If we apply the new LETKF,
several times faster than LEKF, we can further accelerate the process.
3.6 Sensitivity experiments
In this section, we investigate the robustness of our results to the experimental settings.
We apply a difierent observational network and investigate the impact of moisture ob-
servations. Then, we investigate the impact of vertical error correlations and a difierent
type of error covariance localization. Finally, we apply random perturbation addition to
ensemble members. The sensitivity results also suggest possible ways to improve the fllter
performance.
100
3.6.1 Response with a difierent observational network
So far, we applied mostly the sparse observational network, and we saw lattice-like patterns
in the analysis error flelds. Since the observations were available uniformly in space, the
error flelds looked uniform. However, in reality we have more rawinsonde observations
over land and less over the ocean. Thus, we expect non-uniform distribution of error
flelds. In addition, it is not clear how sensitive each data assimilation method is to the
choice of observational network. We expect EnKF to be more robust with respect to the
choice of the observational network than 3DVAR. To see the response with a difierent
observational network, we apply a realistic observational network (Fig.3.11) to 3DVAR
and serial EnSRF as a representative of EnKF methods.
Fig.3.29 shows analysis RMSE of 500hPa height fleld in the cases of 3DVAR (short
dashed line) and serial EnSRF (solid line) with 30 ensemble members and localization
scale of 3.0. Clearly serial EnSRF outperforms 3DVAR. The temporal averages are 45.6
and 11.6 meters for 3DVAR and serial EnSRF, respectively.
Fig.3.30 shows the spatial distribution of analysis RMSE of 500hPa height flelds
of 3DVAR (top panel) and serial EnSRF (bottom panel). Both flgures have a similar
structure with difierent amplitudes. Larger errors appear in data-poor regions. The
southern hemisphere shows larger errors than the northern hemisphere. The largest error,
except the polar area, is observed over south of the Paciflc ocean. The Northern Atlantic
ocean contains the largest error part in the northern hemisphere, except the polar area.
Clearly 3DVAR shows larger errors everywhere, the efiect of the irregular observational
network is larger in 3DVAR. As for the lattice-like pattern, we can see some high frequency
components especially in the tropics in 3DVAR, but we cannot see clear correspondence
to the observational locations. Thus, the error flelds look more realistic; data-rich regions
101
Figure 3.29: Analysis RMSE of 500hPa height fleld in the case of the realistic observational
network using 3DVAR (short dashed line) and serial EnSRF (solid line). The ensemble
size for serial EnSRF is 30, the localization scale is 3.0.
102
show smaller errors, data-poor regions show larger errors.
Figs.3.31 and 3.32 show zonal structure of analysis RMSE of u-wind and temperature
flelds in the cases of 3DVAR and serial EnSRF, respectively. In all flgures, it is clear
that the southern hemisphere has much larger errors than the northern hemisphere. By
comparing 3DVAR and serial EnSRF, both flgures show very similar structure with much
larger amplitudes in 3DVAR. However, there is a structural difierence between 30 and 60
degrees north in the lower to middle atmosphere. Here 3DVAR shows larger errors, serial
EnSRF does not, even if we consider the difierence in the error amplitudes. This is clearer
in the temperature fleld (bottom panel). Although many observations are available in
the northern hemisphere, there are data-poor regions over the Atlantic and Paciflc ocean.
Since 3DVAR is more sensitive to the lack of data, 3DVAR tends to be afiected more over
the data-poor regions. This could explain why 3DVAR shows larger errors in the northern
hemisphere. RMSE of temperature flelds shows large errors over the Eastern Atlantic
especially at lower levels, no other signal such as orographic shape is observed.
3.6.2 Usefulness of moisture observations
Moisture flelds are strongly related to highly nonlinear and complicated physical processes,
moisture flelds tend to show a noisy structure compared to other variables. Since the
moisture flelds represent a complicated and noisy physical processes, it is not clear if the
moisture observations actually provide useful information in data assimilation, especially
in the case of sparse observational networks. In fact, Szunyogh et al. (2005) assumed
moisture was not observed in their LEKF experiments using the NCEP global model. In
spite of these facts, we included moisture observations in the present research, considering
the SPEEDY model has much simpler physical processes than state-of-the-art models. It
103
Figure 3.30: Spatial distribution of analysis RMSE of 500hPa height flelds (in meters) in
the case of the realistic observational network using 3DVAR (top panel) and serial EnSRF
(bottom panel). The ensemble size for serial EnSRF is 30, and the localization scale is
3.0. Dots indicate the observational locations.
104
Figure 3.31: Zonal structure of analysis RMSE of u-wind fleld (top panel) and temperature
fleld (bottom panel) in the case of 3DVAR with the realistic observational network. Shades
and contours show analysis RMSE and mean analysis flelds, respectively.
105
Figure 3.32: The same flgure as Fig.3.31 but in the case of serial EnSRF with 30 ensemble
members and localization scale of 3.0.
106
is important to investigate the impacts of moisture observations, so that we conflrm our
results are not strongly afiected by the noisy moisture observations.
In order to investigate the impact of moisture observations, we perform data assim-
ilation cycle experiments excluding moisture observations and compare with the previous
results including moisture observations. Fig.3.33 shows analysis RMSE using serial En-
SRF with 10 ensemble members with and without moisture observations. For the moisture
variable (speciflc humidity, q), we obtain large improvement by moisture observations as
expected. Importantly, we see positive impact by the moisture observations for other
variables. We also see improvements by moisture observations using LEKF, supporting
moisture observations have positive impact. Thus, it is reasonable to include moisture
observations in our experiments, suggesting our results are not adversely afiected by the
moisture observations.
Note that the SPEEDY model has only very simple physics, whereas the physical
processes in more sophisticated models are much more complicated. As a result, it is
not straightforward to determine in the operational context if moisture observations cause
such big improvement.
3.6.3 Vertical error correlations
So far, vertical error correlations were not considered in any data assimilation method
(3DVAR, serial EnSRF, and LEKF). However, it is known that atmosphere has an impor-
tant vertical structure, especially with respect to synoptic instability such as baroclinic
instability. Thus, we expect that the inclusion of the vertical error correlation might
contribute a certain amount in data assimilation.
To see the impact of vertical correlation, vertical correlation is introduced in LEKF
107
Figure 3.33: Analysis RMSE of the u-wind fleld (a), hight fleld (b), temperature fleld (c)
and speciflc humidity fleld (d) in the cases of serial EnSRF with 10 ensemble members
with and without humidity observations shown by broken and solid lines, respectively.
The thin solid line shows observational error standard deviation wherever applicable.
108
by taking local patches combined with their closest levels. For example, the bottom level
is combined with the second level, thus, the local patch becomes twice larger at the bottom
level, which is the same at the top level. The second level is combined with the bottom
and the third levels, thus, the local patch becomes three times larger at the second level.
Fig.3.34 shows the impact of the vertical correlation in the cases of 10 and 20 en-
semble members. The local patch parameters are flxed to l = 1 and l2 = 0, respectively.
With 10 ensemble members, the vertical correlation has a clear positive impact: 46.82m
RMSE of 500hPa height in average with the vertical correlation (black solid line), while
60.38m without the vertical correlation (black short-dashed line). Thus, the vertical corre-
lation indeed shows a positive impact. With 20 ensemble members, however, it has slight
negative impact in average: 38.35m RMSE of 500hPa height with the vertical correlation
(blue dotted line), while 35.45m without the vertical correlation (blue long-dashed line).
However, the difierence is not statistically signiflcant.
At the other extreme, we can construct the LEKF with full vertical correlations,
i.e. without localization in the vertical. In this case, we combine all the vertical levels
and take a vertical column as analyzed points in a local patch. So far, we needed a three-
dimensional loop for the local patch center to cover the entire three-dimensional grids.
However, now with the full vertical correlation, we need only a two-dimensional horizontal
loop.
Fig.3.35 shows the impact of the full vertical correlation in the cases of 10 and 20
ensemble members. Here, the local patch parameters are chosen to be l = 1 and l2 = 0,
and the line legends are the same as the previous Fig.3.34. With 10 members, the full
vertical correlation shows a clear negative impact. A possible reason is as follows: The
more combined vertical levels, the larger the dimension of a local patch becomes, which is
109
Figure 3.34: Impact of the vertical correlation in the analysis RMSE of 500hPa height fleld
using LEKF with 10 (black lines) and 20 (blue lines) ensemble members. The black solid
line and black short-dashed line show the cases with and without the vertical correlation,
respectively. The blue dotted line and blue long-dashed line show the cases with and
without the vertical correlation, respectively.
110
Figure 3.35: The same flgure as Fig.3.34 but for the full vertical correlation.
why sampling errors in far away levels due to the limited ensemble size get larger. With 20
members, the impact of the full vertical correlation is not clear. It may be possible that the
vertical correlation provides useful information, but sampling errors are also introduced,
even with 20 ensemble members. Thus, the positive and negative efiects are cancelled out
and the net efiect is almost zero.
3.6.4 Error covariance localization
Since the serial EnSRF with 10 ensemble members has outperformed 3DVAR, the ensemble
size of 10 is large enough for EnKF to work in the given conflguration. However, LEKF
with 10 ensemble members has performed worse than 3DVAR. The difierence between
EnSRF and LEKF is in the way of localization. Serial EnSRF localizes error covariance
around observations with Gaussian-like weighting, LEKF localizes around an analyzed
point with Heaviside-like step weighting. Thus, we applied observational error covariance
111
Localization scale ( ) 1 1.0 1/p2 0.5 0.25
500Z RMSE [m] (l = 1) 60.4 48.2 34.5 25.9 42.0
500Z RMSE [m] (l = 2) 75.9 29.5 21.6
Table 3.13: Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble members
with various localization scale of observational error covariance localization. The local
patch parameter is chosen to be l = 1,2.
localization in LEKF.
Table 3.13 shows analysis RMSE of 500hPa height fleld using LEKF with 10 en-
semble members and local patch parameter l = 1 with various localization scales of ob-
servational error covariance localization. We can see a large improvement of LEKF with
the observational error covariance localization in the case of 10 ensemble members. If we
choose localization scale of 0.5, the analysis RMSE (25.9) is smaller than that of 3DVAR
(31.1). When we increase local patch parameter l to 2, we observe better performance
than in the case of local patch parameter l = 1. Thus, by increasing the local patch size,
we include distant observations that still provide useful information.
With 30 ensemble members, LEKF showed 8.44m RMSE of 500hPa height when
local patch parameter l = 4 and 8.67m when l = 3. Both of them are much smaller than
3DVAR. The localization is applied for 30 ensemble members to see the impact of the
localization when the fllter works very well. We chose l = 3 and  = 1.5, then the RMSE
became 7.94m, slightly better than the case without localization.
An important characteristic is observed in the temporal changes of analysis RMSE.
Fig.3.36 shows the temporal sequences of analysis RMSE of 500hPa height flelds with
difierent localization settings. After the initial one-month transient period, the three lines
112
Figure 3.36: Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble mem-
bers. Here, observational error covariance localization is applied. The solid, long-dashed
and short-dashed lines correspond to no localization, with localization scales  = 1.0, and
 = 0.5 respectively.
clearly separate. The order of the performance is already described in Table 3.13. In the
very initial transient period (flrst 10 days or so), the three lines also show clear separation,
but with the opposite order. With large initial errors, distant observations provide useful
information, but the localization reduces their efiect. Thus, the more localized, the smaller
the correction of the initial large error becomes. After the transient period, the sampling
errors in the distant points dominate, therefore, the more localized, the smaller the analysis
errors.
So far, we modifled the localization in LEKF. We can modify the localization of serial
EnSRF using Heaviside-like step weighting function. Houtekamer and Mitchell (1998)
applied this kind of localization around the observations by taking ?radius of in uence?.
113
However, LEKF takes local patch with square shape. Thus, we take a square around an
observation and force zero analysis increment outside of the square so that we simulate
localization weighting function similar to LEKF. The impact of the difierent localization
weighting function turned out to be very small. When we took 1 ? 1 squares around
observations, analysis RMSE of 500hPa height was 14.0m, which is just slightly worse
than 11.7m that was obtained using the more natural Schur product with the length scale
of 1. Here, we used 10 ensemble members.
In EnKF formulations, we applied globally the same localization length scale in
grid spacing. However, as we could see in the background error statistics using the NMC
method, the correlation length scale is not constant (Fig.3.3). We consider the latitudinal
dependence and apply the length scale given by the NMC method in serial EnSRF. Then,
analysis RMSE of 500hPa height was 22.1m, almost twice as large as 11.7m obtained by
the constant localization scale of 1. The large error comes mainly from high latitudes,
suggesting large length scales introduce sampling errors in the estimation of the error
covariance with 10 ensemble members.
3.6.5 Random perturbation addition as ?stochastic seeding?
Corazza et al. (2002) kept the bred vectors (BVs) ?young? by applying stochastic per-
turbations to the BVs used in a 3DVAR data assimilation where the background error
covariance is augmented with BVs in order to account for ?errors of the day?. They found
that the augmentation of the background error covariance with BVs reduced analysis er-
rors by 20%, including stochastic seeding on the BVs at the beginning of the breeding
integration increased the improvement to about 40%. Miyoshi and Kalnay (2005) inves-
tigated the efiect of the ?stochastic seeding? in a breeding cycle on the Lorenz-96 model,
114
Figure 3.37: Analysis RMSE of 500hPa height fleld using LEKF with 10 ensemble mem-
bers. Here, observational error covariance localization is applied. Dashed and solid lines
show the cases with and without random perturbation addition, respectively.
and they found stochastic seeding lets BVs grow secondary instabilities that were ignored
in leading BV but captured in higher order BVs with orthogonalization processes. In
EnKF, the ensemble members do not converge to a single direction, but still, it may be
possible that ensemble members lose signals of instabilities.
Expecting that the random perturbation addition or ?stochastic seeding? has similar
efiects to let EnKF capture possible instabilities, we apply it on LEKF with 10 ensemble
members. Random perturbation consists of random numbers with Gaussian distribution
and 1% of the observational error standard deviation. Fig.3.37 shows analysis RMSE of
500hPa height fleld using LEKF with and without random perturbation addition. With
random perturbation addition, the RMSE is smaller almost everywhere or at least is not
worse than the case without random perturbation addition.
115
3.6.6 Summary
In this section, we conflrmed that our results in the previous section were reasonably
robust with the experimental settings. Using an irregularly distributed realistic obser-
vational network, both 3DVAR and serial EnSRF show larger errors than the regularly
distributed sparse network. 3DVAR is more strongly afiected by the irregular distribution
of observations, the advantage of EnKF is larger. Moisture observations provide useful
information for data assimilation on the SPEEDY model, supporting that we included
moisture observations in our experiments. In LEKF with 10 ensemble members, account-
ing for local vertical correlation improve the results, but full vertical correlation worsens
the fllter. With 20 ensemble members, the efiect of considering vertical correlations is
not signiflcant. These results suggest that the vertical error independence is a reasonable
choice in the SPEEDY model. LEKF shows signiflcant improvements by the observa-
tional error covariance localization, especially with fewer ensemble members. However,
serial EnSRF with 10 ensemble members is not sensitive to the choice of localization
weighting function, suggesting that the shape of localization weighting function does not
have strong in uences. When adding random noise as ?stochastic seeding?, LEKF showed
better results, suggesting possible further improvements in EnKF.
3.7 Characteristics of the analysis and forecast errors
3.7.1 Structures of the analysis increment
Carrying out perfect model experiments, where we know the ?truth?, allows us to ana-
lyze the structure of analysis errors (which are not known in real data assimilation) and
to compare them with forecast errors. It is also important to study the relationship be-
116
tween EnKF estimated errors and observed errors as previously done by Corazza with a
quasi-geostrophic model (Kalnay et al. 2005). We compare forecast errors and analysis
increments for 3DVAR and EnKF.
Fig.3.38 shows the background error flelds (forecast minus truth, shaded) and the
analysis increment flelds (forecast minus analysis, contour) at an arbitrary time in the case
of 3DVAR (top panel) and serial EnSRF (as representative of EnKF) with 30 ensemble
members (bottom panel). Here, the localization scale of serial EnSRF is chosen as 3.0,
which is an optimal. The results we obtain with a primitive-equation model are very
similar to those obtained with the quasi-geostrophic model. The analysis increment fleld
of EnKF captures  ow-dependent structures of the background error fleld, e.g. positive
pressure increment at the lower left corner of the bottom panel. In contrast, 3DVAR
cannot capture the  ow-dependent structures at all. 3DVAR analysis increments look like
circles around observations, and only in some wide error regions with large amplitudes,
increments by neighbor observations look connected, e.g. large positive and negative
increments in higher latitudes. The observations are so sparse that basically 3DVAR
cannot capture the structure between observations.
Importantly, we see a striped pattern in the background error fleld of 3DVAR,
whereas there is no such pattern in the case of EnKF. This is also seen in the zonal RMSE
structure shown in Fig.3.39, where 3DVAR shows lower RMSE every 4 point correspond-
ing to observational location. Although correlations between the analysis increment and
background error are small (around 0.3) both in 3DVAR and EnKF, EnKF shows a clear
advantage because of better background flelds.
117
Figure 3.38: Background error fleld (shaded) and analysis increment fleld (contour) of
surface pressure (ps) at an arbitrary time in the case of 3DVAR (upper panel) and serial
EnSRF with 30 ensemble members and the localization scale of 3.0 (bottom panel).
118
Figure 3.39: Analysis zonal RMSE of 500hPa u-wind flelds in the cases of 3DVAR, serial
EnSRF with 30 ensemble members and the localization scale of 3, and LEKF with 30
ensemble members, the local patch parameter l = 3 and covariance localization scale
 = 1.5.
119
3.7.2 Analysis error and ensemble perturbation growth
5-day forecasts were computed using serial EnSRF and 3DVAR data assimilation cycles.
Here, 10 ensemble members are used in serial EnSRF with the localization length scale
of 1. Fig.3.40 shows perturbation growths of analysis error and ensemble perturbations.
Here, the energy norm (Houtekamer et al. 2005; Ehrendorfer and Errico 1995) is used to
measure the perturbation amplitude:
E = 12
?
u2 +v2 + cpT
r
T2 +RdTr
 p
s
p0
?2!
(3.15)
where E is energy, and constant-pressure heat capacity cp = 1005.7Jkg?1K?1, dry-air gas
constatnt Rd = 287Jkg?1K?1, reference temperature Tr = 270K, and reference pressure
p0 = 1000hPa are constants. E is averaged over three-dimensional grids.
Fig.3.40 shows RMS values of forecast/analysis errors and ensemble perturbations.
?RMSE? shows forecast/analysis RMSE, indicating the growth of the initial analysis er-
rors. ?SPREAD? shows forecast/analysis ensemble perturbation spreads, indicating the
growth of ensemble perturbations. The perturbation amplitudes are averaged in time
using whole February (112 samples) so that the flgure shows one-month averages.
It is noted that 3DVAR analysis RMSE is almost equal to 5-day forecast RMSE
of serial EnSRF. Thus, serial EnSRF provides much better analysis and forecast than
3DVAR in the perfect model case.
Both 3DVAR and serial EnSRF analysis errors grow as expected, and the ensemble
perturbations of serial EnSRF also grow. The question of whether ensemble perturbations
generated by EnKF data assimilation grow or not is not trivial. Houtekamer et al. (2003;
2005) indicated the ensemble perturbations in their EnKF scheme decay in the flrst few
days. However, in their EnKF scheme, they added a model error covariance term to each
120
 0
 5
 10
 15
 20
 25
 30
 35
 40
 0  24  48  72  96  120
ENERGY [m
2
/s
2
]
TIME [HOUR]
SPEEDY-3DVAR/EnSRF 5-DAY FORECAST
RMSE(3DVAR)
SPREAD(EnSRF)
RMSE(EnSRF)
Figure 3.40: The growth of the initial analysis errors and ensemble perturbations
in the cases of 3DVAR and EnKF. Thin solid line shows the growth of the initial
3DVAR analysis errors (RMSE(3DVAR)). Thick solid line and thick broken line show the
growth of the initial EnKF analysis errors (RMSE(EnSRF)) and ensemble perturbations
(SPREAD(EnSRF)), respectively.
121
ensemble. The perturbations due to the model error covariance term are obtained by the
statistical average and are unrelated to the daily  ow. Thus, we expect these perturba-
tions are not dynamically fast growing perturbations, and because of their relatively large
amplitude, they dominate the initial error decay in Houtekamer?s results. Importantly in
our results, the curves of ensemble perturbation and analysis error growth of EnKF are
almost parallel. Thus, ensemble perturbation and analysis error growths show very similar
growth rates (doubling in about three days), which implies both perturbation flelds are
dominated by similar dynamically growing errors.
The ensemble spread is larger than the analysis error, which may be because of the
choice of covariance in ation parameter, (4% multiplicative in ation). This situation is
opposite to what Houtekamer et al. (2005) indicated in their perfect model experiments
(Fig.4 of Houtekamer et al. 2005). They observed that the ensemble spread is smaller
than the analysis error. In the present conflguration, it seems 4% in ation is too large,
which is why the ensemble spread is larger. To conflrm this anticipation, we applied online
estimation of the in ation parameter to satisfy eq.(2.82) (cf. Section 2.3.6). We observed
ensemble spread 2.07 and analysis error 2.425 in the energy norm (originally the spread
was 4.85 and analysis error was 2.25). Thus, if we choose smaller covariance in ation, the
analysis ensemble spread becomes smaller, the relative size between the analysis ensemble
spread and analysis error can be adjusted through the covariance in ation parameter.
However, although the ensemble spread was decreased to the analysis errors by changing
the in ation parameter, the analysis error got slightly larger.
122
3.7.3 Error flelds and ensemble spread
Since EnKF works as expected, ensemble members are expected to capture the error
structures. Thus, we plot the error flelds and ensemble spread in the case of serial EnSRF
with 30 ensemble members and the Schur product length scale of 3.0. The top panel
of Fig.3.41 shows the analysis error fleld at an arbitrary time, and the bottom panel
shows the forecast error fleld at the same time. The analysis and forecast error flelds are
deflned as analysis/forecast minus the true fleld. The shaded fleld in Fig.3.41 shows the
analysis/forecast error fleld, the contour line shows the analysis/forecast ensemble spread.
The ensemble spread represents the analysis/forecast error standard deviation captured
by EnKF. Some areas (indicated by red rectangles) show similar structures between the
error fleld and the ensemble spread. However, the correspondence is not very good, the
pattern correlations between error flelds and ensemble spread are 0.25 in the analysis and
0.18 in the 6-hour forecast. By contrast, the pattern correlation between analysis and
forecast ensemble spreads are very high (0.96). Similarly, the pattern correlation between
analysis and forecast error flelds are also high (0.74).
3.7.4 Forecast errors, analysis errors, and bred vectors
In order to see the shape of the error flelds, we plot the 2-day forecast error and analysis
error flelds valid at the same time. Fig.3.42 shows 2-day forecast error fleld (shaded) and
analysis error fleld (contour) of u-wind at the 4th level ( = 0.51) at an arbitrary time in
the case of serial EnSRF with 10 ensemble members and the localization scale of 1. It is
clear that the two flelds have extremely similar shapes, conflrmed by their high pattern
correlation of 0.76. 2-day forecast error fleld is regarded as a dynamically growing error,
the analysis fleld contains growing errors. Even with EnKF that considers dynamical error
123
Figure 3.41: The error flelds (ananlysis/forecast minus truth; shaded) and the ensemble
spread of 500hPa height fleld (contour) at an arbitrary time. The top and bottom panels
show analysis and background (6-hour forecast), respectively, both at the same time. Some
parts have similar structures indicated by rectangles.
124
Figure 3.42: 2-day forecast error fleld (shaded) and analysis error fleld (contour) of u-wind
at the 4th vertical level ( = 0.51) at an arbitrary time in the case of serial EnSRF with
10 ensemble members and the localization scale of 1. The correlation between the two
flelds are 0.76.
structure, the analysis flelds still contain growing errors. Fig.3.43 shows the same flgure
in the case of 30 ensemble members, which still shows high correlation (0.79). Fig.3.44
shows the same flgure in the case of 3DVAR, which also shows high correlation (0.74).
The three flgures show the same time with the same truth, so it is expected the
dynamically growing error has similar shape. However, all the three flgures show difierent
shapes, which implies there are many possible dynamical ?errors of the day? rather than
a dominance by bred vectors.
One can observe high correlation of the two flelds not only at this moment for the
speciflc variable, but also at every time for all variables. Fig.3.45 shows temporal mean
of the correlation coe?cients for four model variables (u, v, T, q), all of which show high
125
Figure 3.43: The same as Fig.3.42, but in the case of serial EnSRF with 30 ensemble
members and the localization scale of 3. The correlation between the two flelds are 0.79.
Figure 3.44: The same as Fig.3.42, but in the case of 3DVAR. The correlation between
the two flelds are 0.74.
126
Figure 3.45: Temporal mean of correlation coe?cients between 2-day forecast error and
analysis error flelds for four model variables (u, v, T, q) in the case of serial EnSRF with
10 ensemble members and the local patch parameter l = 1.
correlations around 0.75.
In order to further examine the shape of the error flelds, 5-day and 2-day forecast
errors, analysis error, E-dimension, and bred vector (BV) flelds are shown in Fig.3.46.
E-dimension, initially proposed by Patil et al. (2001), is a measure of substantial dimen-
sionality of the space spanned by ensemble members, the deflnition is shown in the next
section. Patil et al. suggested the usefulness of low E-dimensionality. BVs are generated
by a breeding cycle starting from a fleld consisting of random numbers. The time interval
of the breeding cycle is chosen as 6 hours, the rescaling is 3 m2/s2 with the energy norm
(eq.(3.15)), though BVs are robust to the choice of norm. Fig.3.46 (g) shows E-dimension
computed in a 5 ? 5 local patch for a single variable without considering inter-variable
correlations, (h) shows BV. There are some areas showing low E-dimension and large sig-
127
nals in BV such as the north west coast of the United States, the north western corner of
Paciflc ocean, western Siberia, and the southern end of south America. We flnd similar
signals in 5-day forecasts (a) and (b) in these areas. In the 2-day forecasts (c) and (d),
the signal in western Siberia is very clear, but others are not so clear. 2-day forecast may
not be long enough to form growing modes. In the analysis error flelds (e) and (f), we see
a similar signal in the western Siberia with smaller amplitudes. The growing errors shown
in BV and forecast flelds almost disappear in analysis error flelds.
3.7.5 E-dimension and explained variance
Szunyogh et al. (2005) have shown a strong anti-correlation between E-dimension and
explained variance. E-dimension ` is a measure of substantial dimensionality of the space
spanned by ensemble members, deflned as
` = (
Pm
i=1 ?i)2P
m
i=1 ?2i
(3.16)
where m and ? denote the ensemble size and eigenvalues of an m ? m matrix E>E,
cf. eq.(1) of Patil et al. (2001). Here, E is an N ? m matrix consisting of ensemble
members. Explained variance is a measure showing how well the true error is spanned by
the ensemble members. Let G denote a matrix consisting of eigenvectors of the matrix
E>E, explained variance ? is deflned as
? = jjGG
>(xf ?xt)jj2
jj(xf ?xt)jj2 (3.17)
The top and bottom panels of flg.3.47 show zonal averages of the time-mean E-dimension
and explained variance, similar to Figs.5 and 6 of Szunyogh et al. (2005). Here, E-
dimension is computed using a 3 ? 3 local patch for all variables, inter-variable corre-
lations are considered. Similar structure as Szunyogh et al. (2005) is observed: low
128
Figure 3.46: 5-day forecast error ((a) and (b)) and 2-day forecast error ((c) and (d)),
analysis error ((e) and (f)), E-dimension (g), and bred vector (h) of u-wind at the 4th
level all valid at the same time chosen arbitrarily, in the cases of serial EnSRF with 10
members and 3DVAR.
129
E-dimension and high explained variance appear in higher latitudes, high E-dimension
and low explained variance appear in tropics, especially at middle levels and high level.
Fig.3.48 shows the scatter plot between the zonal mean of the time-mean E-dimension
and explained variance, indicating high anti-correlation, similar to Fig.7 of Szunyogh et
al. (2005).
Without the zonal mean, Fig.3.49 shows the time-mean E-dimension (top panel)
and explained variance (bottom panel) at the 4th level ( = 0.51). Fig.3.50 shows the
scatter plot, indicating a high anti-correlation of -0.66. The zonal mean is not essential to
obtain the high anti-correlation between E-dimension and explained variance.
However, if we flx the time, we do not observe the strong anti-correlation. Fig.3.51
shows the same flgure as Fig.3.49 but at a flxed time. Fig.3.52 shows the scatter plot be-
tweenE-dimensionandexplainedvarianceshowninFig.3.51, indicatinglowanti-correlation
of -0.26. There are a few points with low E-dimension and small explained variance, which
explains the weak anti-correlation. If we focus on several regions showing especially low
E-dimension, e.g. east of Australia, north-eastern Indian ocean, east coast of middle South
America, and northern Paciflc, we flnd correspondence of low E-dimension and large ex-
plained variance. In these several parts, EnKF could capture the true error shape and
contribute to reduce the growing errors; we saw growing errors indicated by BV appear
in such areas in the previous section.
3.7.6 Summary
? The analysis increment fleld of EnKF captures  ow-dependent structures of the
background error fleld, 3DVAR cannot capture the  ow-dependent structures at all.
? 3DVAR showed a striped pattern in error flelds, they almost disappeared in EnKF.
130
Figure 3.47: Zonal mean of the E-dimension (top panel) and explained variance (bottom
panel) in the case of serial EnSRF with 10 ensemble members temporally averaged for one
month (112 samples).
131
Figure 3.48: Scatter plot of the zonal mean of the time-mean E-dimension and explained
variance, indicating high anti-correlation.
? Ensemble perturbations generated by serial EnSRF grow with a rate similar to the
growth rate of the analysis error.
? Error and ensemble spread flelds show similar signals in a few areas, but with low
pattern correlations of about 0.2.
? 2-day forecast error and analysis error flelds have very similar shapes, with cor-
relations about 0.75. However, the error flelds are difierent with difierent data
assimilation schemes, suggesting the possibility of multiple ?errors of the day?.
? 5-day forecast error, 2-day forecast error and bred vector flelds show similar signals
in low E-dimensional areas.
? Time-mean flelds of explained variance and E-dimension are anti-correlated with pat-
tern correlation of about -0.66. Without the time-mean, only a few low E-dimension
132
Figure 3.49: The E-dimension (top panel) and explained variance (bottom panel) at the
4th level in the case of serial EnSRF with 10 ensemble members temporally averaged for
one month (112 samples).
133
Figure 3.50: Scatter plot of the time mean E-dimension and explained variance at the 4th
level, indicating high anti-correlation -0.66.
areas show good anti-correlation, with global pattern correlation about -0.26.
3.8 Discussion
In this chapter, three methods (3DVAR, serial EnSRF, LEKF) are implemented and
applied on the SPEEDY model under the perfect model assumption. Our results show
that the EnKF clearly outperforms 3DVAR. One surprising result is that for difierent
methods and number of ensemble members the 2-day forecast ?errors of the day? are very
similar to the analysis errors. However, they are not necessarily similar among difierent
methods, suggesting there are multiple growing ?errors of the day?. In low E-dimensional
regions, however, the errors tend to have similar structures. Bred vectors also show similar
structures in those regions. The initial analysis errors and ensemble perturbations show
very similar growth rates, suggesting the flelds contain growing ?errors of the day?.
134
Figure 3.51: The E-dimension (top panel) and explained variance (bottom panel) at the
4th level in the case of serial EnSRF with 10 ensemble members at an arbitrary time.
135
Figure 3.52: Scatter plot of E-dimension and explained variance at the 4th level at an
arbitrary time, indicating anti-correlation -0.26.
Overall, our results suggest that serial EnSRF outperforms LEKF, but this advan-
tage is substantially reduced when we localize the observational error covariance or increase
the number of ensemble members. A possible reason of the disadvantage of LEKF is that
there are many possible ways to combine overlapped local patches, our implementation of
LEKF may not be optimal. We average the overlapped local analyses in the smaller local
patch denoted by l2 using Heaviside-like step weighting. We see no clear dependencies on
l2. A natural way is to use Gaussian-like weighting according to the distance from the
center point. Further investigation is required to optimize the way to combine overlapped
local patches in the LEKF. It is possible that the formulation of the local patches over
the poles may not be optimal in the LEKF.
In a realistic case with a parallel computer and a large number of observations,
LEKF should be much faster than the serial EnSRF, especially with the new LETKF
136
developed by Hunt (pers. comm.). As a result, in an operational setup, LEKF/LETKF
may be the only computationally feasible choices. However, in the present experiments,
we used a single computer and a small number of observations. For this case, LEKF
has shown no advantage against serial EnSRF; serial EnSRF performs better with less
computational time. Thus, we choose serial EnSRF for the last part of the experiments.
137
Chapter 4
Data assimilation in the presence of model errors
4.1 Introduction
4.1.1 Overview
In Chapter 3, three data assimilation methods (3DVAR, serial EnSRF, and LEKF) are
applied on the SPEEDY model. EnKF clearly outperforms 3DVAR under the perfect
model assumption. In this chapter, we remove the perfect model assumption using the
NCEP/NCAR reanalysis (NNR) flelds (Kalnay 1996) as the nature run. We explore the
following questions in this chapter:
1. What do model error bias and EOFs (empirical orthogonal functions, e.g. Wilks
1995) look like?
2. What are the efiects of model errors on data assimilation?
3. How does model bias estimation (Dee and da Silva 1998) work with EnKF?
We introduce NNR and compute model error statistics by taking difierences between
NNR flelds and the SPEEDY model 6-hour forecasts initialized by NNR flelds. Then, the
efiects of model errors on data assimilation are investigated. We perform data assimilation
cycle experiments using the same data assimilation systems as in Chapter 3, but replacing
the nature run by NNR flelds. We describe theory on model bias estimation following Dee
and da Silva (1998) and apply the estimation method with EnKF. We also apply low-order
bias estimation, assuming the bias fleld is spanned by a limited number of base flelds. As
base flelds, we use the model error bias fleld and EOFs.
138
4.1.2 OSSE methodology
So far, in this research, we followed OSSEs (Observing Systems Simulation Experiments)
as described in sections 1.3 and 3.4. In OSSE methodology, we generate observations using
the nature run generated by the model forecasts. So far, we have used the same model for
both generating observational data and performing data assimilation. This perfect model
experiment is over optimistic partly because such observational data do not contain the
error of representativeness. Alternatively, Atlas (1997) suggested to generate observations
using a more sophisticated and higher resolution model, difierent from the model used in
data assimilation. This considers the error of representativeness as well as the imperfection
of the model. Similarly in this chapter, we use NNR flelds to generate observational data,
the NNR flelds are generated by a relatively recent sophisticated model (1995 NCEP
model) and data assimilation system. The difierence from the Atlas?s approach is that
we consider the real nature. Real observations are assimilated in NNR flelds, NNR flelds
simulate the evolution by the unknown ?true nature? model.
4.1.3 NCEP/NCAR reanalysis (NNR)
The NCEP/NCAR reanalysis (NNR) flelds (Kalnay 1996; Kistler et al. 2001) are a fairly
good estimate of the true atmospheric state for the past 50 years. They are generated by
a relatively recent (1995) sophisticated model and data assimilation system. We assume
the NNR flelds as a proxy for the truth. The validity of this assumption is beyond the
scope of this research. We randomly perturb the NNR flelds to simulate observations with
the same standard deviation (Table 3.2), we verify our analysis flelds against the NNR
flelds.
NNR data are freely available on the Internet. The NNR data on pressure levels
139
have a difierent grid structure (144?73?17) than the SPEEDY model (96?48?7). The
NNR flelds need to be interpolated to obtain the flelds on the SPEEDY model grids. The
horizontal interpolation is a flrst order linear interpolation. No vertical interpolation is
required to obtain pressure level data because the reanalysis data have all seven pressure
levels: 925, 850, 700, 500, 300, 200, and 100hPa. A flrst order linear interpolation is
applied to obtain sigma level data. Using horizontally interpolated ps of the reanalysis
data, ps ? = p is computed to get the target pressure level for the interpolation. The
reanalysis fleld is interpolated for the target level. logp is used as a vertical coordinate
for the linear interpolation.
Since the original NNR flelds have a higher resolution, they contain wavenumber
components higher than those which the SPEEDY model can resolve. Wavenumber com-
ponents higher than 30 cannot be detected in the SPEEDY model since it has spectral
triangular truncation of 30 wavenumbers. The reanalysis flelds are smoothed out by the
spectral transformation of the SPEEDY model.
For the present experiments, we downloaded the NNR flelds from the NCAR?s web-
site and generated the interpolated and smoothed flelds during the period from January
1, 1982 to February 28, 1982. This period coincides with the satellite era that started in
1979. For simplicity, hereafter we call the interpolated and smoothed NNR flelds ?reanal-
ysis flelds?.
4.2 Model error statistics
In this section, model error statistics (bias and EOFs) are computed and shown.
140
4.2.1 Model error biases
In order to estimate the model error against the reanalysis flelds, the reanalysis flelds are
given as initial conditions to the SPEEDY model. Since the reanalysis flelds are evolved
by an unknown model expected to simulate the true atmospheric dynamics, the difierences
between the SPEEDY model 6-hour forecasts and reanalysis flelds valid at the same time
represent the model errors. We assume the difierence is a sample of model errors.
Fig.4.1 shows the difierences of the three prognostic variables (u, v and T) between
the SPEEDY model 6-hour forecasts and reanalysis flelds, temporally averaged for two
months using 235 samples in the period from January 1, 1982 to February 28, 1982.
Temporally averaged difierences are interpreted as the model bias. Large areas are not
shaded by colors, suggesting that the SPEEDY model simulates the real atmosphere fairly
well. The largest model bias of the wind flelds can be seen in polar regions. Some
difierences are caused by the orographic efiect, most clearly seen in the u-wind fleld at
200hPa level, where there is a clear negative bias downstream of the Himalayas. As
for temperature flelds, a strong negative bias is observed at lower levels over Antarctica,
though almost no clear bias can be seen at 500hPa level. Strong biases are observed mostly
over land at lower levels.
Fig.4.2 shows the same flgure as Fig.4.1 but for the variables q (speciflc humidity), z
(height) and ps (surface pressure). q biases are large in the tropics. Since z is a diagnostic
variable obtained from ps and T, it is clear that z biases are strongly afiected by the ps
bias fleld. It seems that the orographic difierence in the reanalysis fleld and the SPEEDY
model afiects the ps model bias. However, it still shows some biases over the ocean, e.g.
positive bias over the western tropical paciflc basin, which disappears in the z bias fleld
in the upper troposphere.
141
Figure 4.1: SPEEDY model biases of u, v [m/s] and T [K] relative to the NCEP/NCAR
reanalysis flelds.
142
Figure 4.2: SPEEDY model biases of q, z [m] and ps [Pa] relative to the NCEP/NCAR
reanalysis flelds. Each vertical level has difierent scales for q flelds, it is omitted. Red
color shows larger values.
143
4.2.2 EOFs
We computed empirical orthogonal functions (EOFs, see for example, Wilks 1995) of
the two-month model error samples. EOFs are deflned as eigenvectors of the covariance
matrix C of the model error anomaly x0(t,s) where t and s denote time and space indices,
respectively:
C = x0(t,s)>x0(t,s)
= VDV> (4.1)
The model error anomaly x0(t,s) is deflned by subtracting the time-mean model bias fleld
?x(s) from time series of model error flelds x(t,s):
x0(t,s) = x(t,s)? ?x(s) (4.2)
V is an N?N matrix composed of eigenvectors, D is an N?N diagonal matrix composed
of eigenvalues, where N denotes the number of the space index. Usually, N is a number
of O(104) or larger, the number of time index m is smaller. As a result, the covariance
matrix is degenerate. Ignoring zero eigenvalues, we have V as an N ?m matrix and D
as an m?m diagonal matrix. We write the singular value decomposition of x0(t,s)
x0(t,s) = UD1/2V> (4.3)
where U is an m?m matrix composed of left singular vectors, D1/2 is an m?m diagonal
matrix composed of singular values, and V is an N?m matrix composed of right singular
vectors. Since singular vectors are orthonormal, U>U = I is satisfled. Using eq.(4.3), we
get EOFs as follows:
144
x0(t,s)>x0(t,s) =
?
UD1/2V>
?>?
UD1/2V>
?
= VD1/2U>UD1/2V>
= VDV> (4.4)
V is equal to that in eq.(4.1), containing EOFs. There is an algorithm to flnd the largest
singular values and corresponding singular vectors, we obtain flrst few EOFs with a little
computation.
Here, we take space index s by combining variables (u, v, T, and q) but separating
levels. In this way, inter-variable correlations are considered but no correlations between
levels are considered. Fig.4.3 shows the flrst EOF (top panel) and the second EOF (bot-
tom panel) of temperature flelds at the bottom level (shades) and surface pressure flelds
(contour). There is a clear signal corresponding to diurnal variation with wave number 1
over longitudes. We also see the phase difierence between the flrst and second EOFs. The
SPEEDY model considers only daily updated forcings and no diurnal variation. Temper-
ature at the bottom level is strongly afiected by insolation, showing strong diurnal signals
over land. Other variables at other levels do not show this kind of clear structures. The
explained variances are 20%, 15%, 12%, 10% in the flrst four EOFs. Fig.4.4 shows time
series of expansion coe?cients of the flrst four EOFs for a week. We see a clear signal of
the diurnal cycle for the flrst and second EOFs. The second EOF shows diurnal  uctu-
ation with a slightly smaller amplitude and a difierent phase than the flrst EOF. Fig.4.5
shows time series of expansion coe?cients of the flrst EOF (top panel) and the third and
fourth EOFs (bottom panel) in a longer time range. Diurnal  uctuation is dominant in
the flrst EOF, but we also see longer time  uctuation with a small amplitude. We see
145
a very similar time series for the second EOF. In the third and fourth EOFs, there are
longer time  uctuations spanning positive and negative values with periods of about 1-2
weeks. Since the training period was only 50 days, it is clearly insu?cient to estimate in
a robust way variable errors with a time scale not much smaller than the training period.
Although the explained variances do not show a clear jump after the second EOF, the
flelds themselves and expansion coe?cients indicate that with this short training period,
only the flrst two EOFs (representing the error due to the absence of diurnal cycle) can
be estimated in a statistically signiflcant fashion.
4.3 Efiects of model errors on data assimilation
While the previous section showed bias and EOFs of the model errors, the efiects of sys-
tematic model errors on data assimilation are not clear. To see the efiects, we perform the
same data assimilation experiments as the perfect model experiments using the reanalysis
flelds as truth. In other words, the same data assimilation systems are applied, the only
difierence is the nature run is replaced by the reanalysis flelds. In this chapter, only serial
EnSRF is used as representative of EnKF as described in Chapter 3, we will refer to it
as ?EnKF? hereafter. The parameter values are chosen as optimal, shown in Table 3.9,
Chapter 3.
Fig.4.6 shows the analysis RMSE of 500hPa height fleld. Even EnKF with 30 en-
semble members, the best among the three, shows about 60 meters of analysis RMSE.
3DVAR shows the largest temporal  uctuation and largest mean errors. However, the
performance of EnKF is more than 10 times worse than in the perfect model case, and the
ratio of 3DVAR to EnKF is clearly smaller than in the perfect model case. This means
that EnKF is more strongly afiected by the model errors.
146
Figure 4.3: EOFs of two-month model error samples from January 1, 1982 to February
28. The top and bottom panels show temperature flelds at the bottom level (shades) and
surface pressure flelds (contour) of the flrst and second EOFs, respectively.
147
1 2 3 4 5 6 7 8?200
?150
?100
?50
0
50
100
150
TIME[DAY]
EXPANSION COEFFICIENTS
EOF1EOF2
EOF3EOF4
Figure 4.4: Time series of expansion coe?cients of the flrst to fourth EOFs for a week.
Solid, dashed, dotted, dash-dotted lines show flrst, second, third, and fourth EOFs, re-
spectively.
148
5 10 15 20 25 30 35 40 45 50?200
?150
?100
?50
0
50
100
150
TIME[DAYS]
EXPANSION COEFFICIENTS
EOF1
5 10 15 20 25 30 35 40 45 50?200
?150
?100
?50
0
50
100
150
TIME[DAYS]
EXPANSION COEFFICIENTS
EOF3EOF4
Figure 4.5: Time series of expansion coe?cients of the flrst EOF (top panel) and the third
(solid line) and fourth (dotted line) EOFs (bottom panel) for 50 days.
149
Figure 4.6: Analysis RMSE of 500hPa height fleld of 3DVAR (short dashed line) and
EnKF with 10 ensemble members (solid line) and 30 ensemble members (long dashed
line) in the case of the regular sparse observational network (Fig.3.10) in the presence of
model errors.
150
Figure 4.7: Analysis RMSE of 500hPa height fleld of 3DVAR and EnKF with 30 ensemble
members in the case of the sparse and realistic observational networks in the presence of
model errors. Blue lines show the cases of realistic observational network, black lines show
the cases of sparse observational network. As for blue lines, short dashed line and long
dashed line show 3DVAR and EnKF, respectively. As for black lines, solid line and dotted
line show 3DVAR and EnKF, respectively.
When the realistic observational network (Fig.3.11) is applied, the advantage of
EnKF with respect to 3DVAR becomes larger. Fig.4.7 shows the analysis RMSE of 500hPa
height fleld in the cases of both the sparse and realistic observational networks. Clearly,
the difierence between 3DVAR and EnKF is larger in the case of the realistic observational
network. 3DVAR is strongly afiected by the difierent observational network, the temporal
 uctuation is larger with the realistic network. By contrast, EnKF is not sensitive to the
choice of the observational network.
Fig.4.8 shows time-mean RMSEs of other variables at other levels. Importantly,
151
3DVAR and EnKF with 10 ensemble members show similar results. Signiflcant difierence
is only seen in height and temperature flelds at a few levels; EnKF shows advantage in
height flelds at middle and upper levels, 3DVAR shows advantage in temperature flelds
at 500, 200, and 100hPa levels. The fact that EnKF is less sensitive to the choice of
observational network than 3DVAR is generally valid for other variables at other levels
except temperature and speciflc humidity, and height flelds at upper levels. There is
almost no difierence among the flve data assimilation methods in speciflc humidity flelds.
4.4 Model bias estimation
The previous section showed that model errors afiect data assimilation, especially in EnKF.
The advantage of EnKF to 3DVAR became smaller than the perfect model case. In this
section, model bias estimation (Dee and da Silva 1998) is introduced and applied to
improve the fllter performance in the presence of model errors.
4.4.1 Theory
We follow the same approach as Dee and da Silva (1998) for bias estimation within KF.
In Section 2.2, we assumed a perfect model. Model imperfections introduce additional
suboptimality to the data assimilation. Thus, we need to modify eq.(2.18) to include a
model bias:
?xfi+1 = M(?xai )?bai (4.5)
Here, ?x denotes an unbiased state. Since the model M has a bias b, it is necessary to
subtract the bias to get the unbiased forecast. The bias also evolves in time:
bfi = Mb(bai?1) (4.6)
152
Figure 4.8: Analysis RMSE at the all pressure levels temporally averaged for one month af-
ter the initial one-month spin-up period (112 samples) in the cases of 3DVAR (short-dashed
lines), EnKF with 10 ensemble members (solid line), EnKF with 30 ensemble members
(long-dashed lines). Black and blue lines show the cases with the regular sparse observa-
tional network and the realistic observational network, respectively. The four panels (a),
(b), (c) and (d) correspond to u-wind, height, temperature (T) and speciflc humidity (q)
flelds, respectively. The observational error standard deviations are shown as thin solid
lines wherever applicable.
153
where Mb denotes the time evolving model of the bias. Now, we have two separate
variables: the model state variable x and the bias variable b.
For the unbiased variable ?x, we can apply the same flltering equations as in Section
2.2, thus,
?xai = ?xfi +K(yo ?H?xfi ) (4.7)
where the optimal weight K is given as
K = PfH>(HPfH> +R)?1 (4.8)
cf. eq.(2.15). Here, the error covariance matrix Pf is deflned for the unbiased variable as
Pf =
D
(?xf ?xt)(?xf ?xt)>
E
(4.9)
The covariance analysis is given by
Pa = (I?KH)Pf (4.10)
cf. eq.(2.17). With the ensemble formulation, the covariance forecast is derived from the
ensemble forecast, and the ensemble update is given by the SRF (square root flltering) to
satisfy eq.(4.10).
Similarly for the bias, we deflne a weighted mean of the bias forecast and the bias
observation:
bai = bfi ?Kb(yo ?Hxfi +Hbfi ) (4.11)
Similar to eqs.(2.2) and (2.3), we deflne the difierences from the true bias bt:
?ba = ba ?bt (4.12)
?bf = bf ?bt (4.13)
154
Subtracting bt from both sides of eq.(4.11), we get
?ba = ?bf ?Kb
?
?yo ?H?xf +H(?bf +bt)
?
= ?bf ?Kb
?
?yo ?H(?xf ?bt)+H?bf
?
= ?bf ?Kb
?
?yo ?H??xf +H?bf
?
(4.14)
where ?xf is deflned as the same as eq.(2.3) in Section 2.2, and we assumed unbiased
forecast error:
??xf = xf ?xt ?bt (4.15)
= ?xf ?xt (4.16)
= ?xf ?bt (4.17)
Bias analysis error covariance Pab is written
Pab =
D
?ba(?ba)>
E
=
??
?bf +Kb(?yo ?H??xf +H?bf)
??
?bf +Kb(?yo ?H??xf +H?bf)
?> 
=
??
(I?KbH)?bf ?Kb?yo +KbH??xf
??
(I?KbH)?bf ?Kb?yo +KbH??xf
?> 
= (I?KbH)Pfb(I?KbH)> +KbRK>b +KbHPfH>K>b (4.18)
where bias forecast error covariance is deflned as
Pfb =
D
?bf(?bf)>
E
(4.19)
and all the cross covariances are assumed to be zero:
D
??xf(?yo)>
E
= 0 (4.20)
D
?bf(?yo)>
E
= 0 (4.21)
D
??xf(?bf)>
E
= 0 (4.22)
155
In order for the weighted mean eq.(4.11) to be optimal, we minimize the total bias analysis
error variance, i.e. the trace of Pab. Thus,
?
?Kb (trace(P
a
b)) = 0 (4.23)
is satisfled (cf. eq.(2.11)). Using the formula (2.12) and (2.13), we get the optimal weight
Kb = PfbH>(HPfbH> +HPfH> +R)?1 (4.24)
cf. eqs.(40) and (53) in Dee and da Silva (1998). Substituting the optimal weight (4.24)
into eq.(4.18), we can get
Pab = (I?KbH)Pfb (4.25)
cf. eq.(54) in Dee and da Silva (1998).
The dimension N of the bias variable b is the same as that of the model variable x,
and the bias error covariance Pb is an N?N matrix. To avoid the explicit computation of
the large matrix, we need a low dimensional assumption to apply the ensemble formulation.
We expand all the covariances with a limited number of ensemble members in a similar
way as in Section 2.2. Thus, we use a bias ensemble to represent Pb.
Since we do not have good information about the bias forecast model Mb, persistence
is assumed as the bias forecast model:
Mb = I (4.26)
which assumes the model bias should be intrinsically constant in time. Then, eq.(4.6) is
just
bfi = bai?1 (4.27)
In this case, the bias has no dynamics, and the low dimensional assumption may not be
reasonable.
156
Alternatively, we can choose a limited number, say k, of base flelds to span the bias
b (Dee and da Silva 1998; Danforth, Kalnay, pers. comm.):
b = Tfl (4.28)
Here, T is an N ?k matrix where each column denotes base flelds, and fl denotes a k-
dimensional vector. In this way, fl is the variable to be estimated, T is flxed. For the best
estimate of fl, we substitute eq.(4.28) into the previous equations derived for b. Here, the
error covariance of b is substituted as
Pb =
D
?b(?b)>
E
=
D
T?fl(T?fl)>
E
= T
D
?fl(?fl)>
E
T>
= TPflT> (4.29)
Thus, substituting eqs.(4.28) and (4.29), we construct the Kalman flltering equations for
fl as follows:
flfi = flai?1 (4.30)
Pffl = (1+?)Pafl (4.31)
flai = flfi ?Kfl(yo ?Hxfi +HTflfi ) (4.32)
Kfl = PfflT>H>(HTPfflT>H> +HPfH> +R)?1 (4.33)
Pafl = (I?KflHT)Pffl (4.34)
where ? denotes a covariance in ation parameter for Pfl. If the dimension of fl is small
enough, we can compute the above equation explicitly.
In summary, we introduced two bias estimation schemes:
157
1. EnKF for the full dimensional bias variable (Dee and da Silva 1998)
2. KF for the low dimensional bias variable (suggested by Kalnay and by the results of
Danforth et al. 2005)
4.4.2 Numerical experiments
First, we applied the full dimensional bias variable estimate using the ensemble formulation
of Kalman flltering. We use the serial EnSRF system with minimal modiflcations. The
difierence is in the gain computation. Since we need to consider not only Pfl but also Pf
in the gain computation, we cannot apply exactly the same system for the bias variable.
However, because we have the forecast ensemble as well as the bias ensemble in memory,
the modiflcation is not complicated.
The initial bias ensemble is generated using Gaussian random numbers. A random
perturbation fleld is constructed by putting a Gaussian random number at each grid point
and smoothing high frequency components using a Lanczos fllter. The amplitudes are
normalized using observational error size, so that each random variable has an appropriate
order of magnitude. We used 10 ensemble members, where 5 random flelds are generated
and their plus-minus pairs constitute 10 members. Thus, initially the bias is assumed to
be zero, i.e. ensemble mean is zero.
When we applied Dee and da Silva?s full dimensional bias estimation, the fllter was
unstable and diverged rapidly. In order for the fllter to be more stable in the initial
transient period, the bias correction eq.(4.5) is modifled as
?xfi+1 = M(?xai )?fibai (4.35)
where fi is a scalar parameter that weights the bias correction. Even with fi = 0.01, one
percent bias correction, we still observed unstable flltering. Fig.4.9 shows the analysis
158
Figure 4.9: Analysis RMSE of 500hPa height fleld of EnKF with 10 ensemble members in
the case of the sparse observational networks in the presence of model errors. The dashed
and solid lines show the cases with and without Dee and da Silva?s full dimensional bias
estimation, respectively.
RMSE of 500hPa height fleld with/without bias estimation. Here, we used fi = 0.01. The
time period flnishes before February, because the fllter diverges just after the flnal time of
the flgure even with fi = 0.01. In the flrst week, the bias estimate makes analyses better,
but eventually it gets worse.
Possible reasons are as follows:
? Too small ensemble size
? Too few observational stations
As mentioned before, the bias forecast is assumed to have no dynamics. Thus, the bias fleld
cannot have dynamically growing modes. Therefore, the low-dimensional assumption for
159
ensemble formulation may not be a good assumption. Even with localization, 10 ensemble
members may be too few to capture all the possible directions. As for the number of
observations, larger number of observations helps the fllter converge faster. With fewer
observations and without dynamics, the initial spin-up period of the flltering becomes
longer.
Since the full dimensional estimation made the fllter unstable, we applied low-
dimensional bias variable estimation. We used the mean bias flelds shown in Figs.4.1
and 4.2 to span the model bias. Although variables are combined, each level is separated.
For comparison, we show the case with constant bias subtraction, i.e. the bias variable fl
is flxed to 1. In addition, we also show the cases of 3DVAR with and without constant bias
subtraction. Fig.4.10 shows analysis RMSE of 500hPa height fleld. Thin red line shows
6-hour forecast errors initiated from the reanalysis flelds. These forecasts provide back-
ground errors purely caused by model errors. Clearly, EnKF shows smaller analysis RMSE
with bias estimation. There is no signiflcant difierence between low order bias estimation
(green short-dashed line) and constant bias subtraction (blue dotted line). 3DVAR (blue
lines) shows slightly better performance with constant bias subtraction, but the improve-
ment is smaller than with EnKF. Fig.4.11 shows the values of all the seven components
of the bias variable. The seven components correspond to vertical levels. After the initial
spin-up, all values are close to one with small temporal variations which explains why the
constant bias subtraction shows similar performance.
Fig.4.12 shows the spatial distribution of analysis RMSE of 500hPa height flelds
with and without bias estimation. Clearly, the amplitude is reduced more in the case with
bias estimation (bottom panel) than the case without bias estimation (top panel). The
large difierence can be seen for example over the Andes of South America, Himalayas,
160
Figure 4.10: Analysis RMSE of 500hPa height fleld for EnKF with low order bias estima-
tion (green short-dashed line), EnKF with constant bias subtraction (black dotted line),
3DVAR with constant bias subtraction (blue dotted line). Solid lines show the cases with-
out bias estimation in the cases of EnKF (black line) and 3DVAR (blue line), the same as
solid and long-dashed lines of Fig.4.6. Thin red line shows 6-hour forecast errors initiated
from the reanalysis flelds.
161
-0.5
 0
 0.5
 1
 1.5
 40  80  120  160  200
TIME [DA CYCLES]
SPEEDY-EnSRF BIAS VARIABLES
1
2
3
4
5
6
7
Figure 4.11: Values of all the seven components of the bias variable when the mean
bias fleld is used as a basis to span the bias where each vertical level is assumed to be
independent. The seven lines correspond to vertical levels.
162
and southern Persian Gulf, where the top panel shows areas with large errors while the
bottom panel does not. In such areas, large model error biases are evident in Fig.4.2.
To see what caused the large error spots, the analysis error bias flelds are shown in
Fig.4.13. The case without bias estimation (top panel) shows much larger amplitudes than
the case with bias estimation (bottom panel), large error spots exist in similar locations.
If we subtract the error bias from the error flelds and take the RMS, i.e. analysis error
standard deviation, the difierence between the cases with and without bias estimation
becomes small. Fig.4.14 shows that the difierence between the top and bottom panels are
smaller than the previous flgures. Here, the top and bottom panels show the cases without
and with bias estimation, respectively. Importantly, the large error spots disappear in the
top panel. Thus, we deduce that the large error spots are caused by the persistent model
bias, which is greatly reduced by the bias estimation.
Fig.4.15 shows analysis RMSE of other variables at all vertical levels. Within EnKF,
the bias estimation generates better analysis states for all variables at all levels. The dif-
ference is 30-40 meters in height flelds, 2-3 m/s in wind flelds, 1-3 K in temperature
flelds, all of which are fairly large improvements. There is almost no difierence between
the cases with low order bias estimation (green short-dashed line) and with constant bias
subtraction (black dotted line), only speciflc humidity RMSE shows signiflcant difierence.
3DVAR and EnKF show similar performance without bias subtraction (solid lines). How-
ever, with bias subtraction, EnKF shows an advantage especially in the wind and height
RMSEs. In the temperature RMSE, EnKF shows an advantage of less than 0.5 K. In
the speciflc humidity, the situation is opposite, 3DVAR shows smaller RMSE than EnKF.
With low order bias estimation, EnKF shows almost identical performance as 3DVAR
with constant bias subtraction. Concentrating on dynamics, we obtain much better per-
163
Figure 4.12: Spatial distribution of analysis RMSE of 500hPa height flelds (shades in me-
ters) in the case of imperfect model experiments using EnKF with 10 ensemble members.
Contours show mean analysis flelds (in meters). Top and bottom panels show the cases
without and with bias estimation, respectively. Dots indicate the observational locations.
164
Figure 4.13: Analysis error bias flelds of 500hPa height (shades in meters) in the case
of imperfect model experiments using EnKF with 10 ensemble members. Contours show
mean analysis flelds (in meters). Top and bottom panels show the cases without and with
bias estimation, respectively. Dots indicate the observational locations.
165
Figure 4.14: Analysis error standard deviation flelds of 500hPa height (in meters) in the
case of imperfect model experiments using EnKF with 10 ensemble members. Contours
show mean analysis flelds (in meters). Top and bottom panels show the cases without and
with bias estimation, respectively. Dots indicate the observational locations.
166
formance using EnKF with bias estimation and constant bias subtraction in wind and
height flelds. However, for variables strongly related to physical processes, the advantage
of EnKF to 3DVAR is less signiflcant. There is signiflcant improvement of humidity flelds
with low order bias estimation compared to constant bias subtraction, similar to 3DVAR
with constant bias subtraction.
Fig.4.16 shows the zonal structure of analysis RMSE of u-wind and temperature
flelds in the case without bias estimation. Fig.4.17 shows the same flgure but in the case
with bias estimation. Again, the amplitude is larger in the case without bias estimation.
Both cases have similar structures in the zonal mean errors. u-wind error flelds show a
clear striped pattern according to the observational locations. We see a large reduction
of temperature errors in the lower altitudes in the northern hemisphere in the case with
bias estimation.
We include the flrst few EOFs as bases to span the bias variables. Table 4.1 shows
analysis RMSE of 500 hPa height fleld with various numbers of EOFs. The mean bias
fleld and flrst two EOFs give the best performance, in agreement with the argument in
Section 4.2.2 the higher order EOFs do not contain statistically signiflcant information
due to the shortness of the training period. Since the diurnal cycle components and the
bias are well deflned, even with a short training period, their online correction leads to
signiflcant improvements. The addition of EOFs beyond the diurnal cycle (not shown)
makes the results slightly worse.
4.5 Summary and discussion
The bias estimation method proposed by Dee and da Silva (1998) did not succeed for full
dimensionality. Thus, we applied a low-dimensional expansion of bias flelds, as suggested
167
Figure 4.15: Analysis RMSE at the all pressure levels temporally averaged for one month
after the initial one-month spin-up period (112 samples) in the cases of EnKF with low
order bias estimation (green short-dashed line), EnKF with constant bias subtraction
(black dotted line), 3DVAR with constant bias subtraction (blue dotted line). Solid lines
show the cases without bias estimation in the cases of EnKF (black line) and 3DVAR
(blue line), the same as the black solid line and the black short-dashed line of Fig.4.8.
The four panels (a), (b), (c) and (d) correspond to u-wind, height, temperature (T) and
speciflc humidity (q) flelds, respectively. The observational error standard deviations are
shown as thin solid lines wherever applicable.
168
Figure 4.16: Zonal structure of analysis RMSE of u-wind fleld (top panel in m/s) and
temperature fleld (bottom panel in K) in the case of imperfect model experiments using
EnKF with 10 ensemble members where no bias estimation is applied. Shades and contours
show analysis RMSE and mean analysis flelds, respectively.
169
Figure 4.17: The same flgure as Fig.4.16 but with bias estimation.
170
Bias correction Uncorrected Mean bias EOF1 EOF2
500Z RMSE [m] 66.6 42.8 40.6 39.2
Table 4.1: Analysis RMSE of 500hPa height fleld using EnKF with 10 ensemble members
and low order bias estimation using various numbers of bias basic flelds. ?Uncorrected?
denotes the case without bias correction. ?Mean bias?, ?EOF1?, and ?EOF2? denotes the
cases with low-order bias correction using the mean bias fleld, using the mean bias fleld
and flrst EOF, and using the mean bias fleld and flrst and second EOFs, respectively.
by Kalnay and by the results of Danforth et al. (2005). The theory of low-dimensional
expansion was also described in Dee and da Silva (1998).
We computed statistics using model errors sampled by the difierence between the
reanalysis flelds and 6-hour forecasts started from the reanalysis flelds. We computed the
mean bias and EOFs from the model error samples, and used them as bases of model
bias. Then, we found the low order bias estimation has a positive impact even with a
single statistical bias fleld. Constant bias subtraction has also shown similar performance,
although it resulted in larger RMSE in speciflc humidity flelds. Although 3DVAR and
EnKF have similar performance without bias estimation or constant bias subtraction,
constant bias subtraction showed larger positive impact within EnKF than within 3DVAR,
except in speciflc humidity. When we applied low order bias estimation within EnKF,
EnKF outperformed 3DVAR everywhere and for all flelds. It is important to note that
we used the dependent sample for the statistics, suggesting we may have overestimated
the performance. It is expected that if we sample with other terms or years to obtain the
basic bias flelds, our results would be worse than what we show in the present research.
However, other experiments have shown similar basic bias flelds even using sampling of
171
other years (Danforth et al. 2005).
Another problem is that the SPEEDY model has much larger errors than sophis-
ticated operational models. We used the NNR flelds as a best estimate of the true at-
mosphere and took the difierence from the SPEEDY model forecast to obtain samples of
model errors. In practice, we use a state-of-the-art model both for estimating the true
nature (i.e. reanalysis) and for forecasts. If we use the same model for both reanalyses and
forecasts, the difierences between the forecasts and analyses by the same model, i.e. the
analysis increments, will be samples of model errors. The model bias fleld obtained from
the analysis increments was very difierent from the one obtained from the NNR flelds.
Using the analysis increments as model error samples, we obtained worse results (not
shown). Thus, obtaining good samples of model errors may be a big challenge in practice,
requiring longer training periods, and/or a multi-model approach; we may obtain a best
estimate of the true nature by an optimal combination with separate analysis flelds such as
NNR and ECMWF (European Centre on Medium-range Weather Forecasting) reanalysis
flelds.
172
Chapter 5
Summary and future directions
5.1 Summary
Our goal was to develop a path towards the operational application of ensemble Kalman
flltering (EnKF), which is expected to produce substantially improved weather/climate
information. Several approaches to EnKF for atmospheric systems have been suggested,
but were not previously directly compared. In addition, there has been little research
published on the sensitivity of EnKF to the imperfections of forecast models. The present
research aimed to explore two basic questions:
1. What are the relative advantages and disadvantages of the two most promising EnKF
methods?
2. How large are the efiects of model errors on data assimilation, and can they be
reduced by model bias correction?
In Chapter 2 we described a theoretical review of EnKF, followed by the FORTRAN
development and testing on the Lorenz-96 model (Lorenz 1996; Lorenz and Emanuel 1998)
of two EnKF methods: 1) a serial ensemble square root fllter (serial EnSRF, Whitaker and
Hamill 2002) that assimilates observations sequentially and localizes the background error
covariance multiplying it with a Gaussian-like function; and 2) a local EnKF (LEKF, Ott
et al. 2002; 2004) that assimilates all the observations within a local volume surrounding
a grid point simultaneously. We reproduced the results obtained by Whitaker and Hamill
(2002) and Ott et al (2004). With observational error covariance localization (multiplying
173
the observational error covariance by the inverse of the Gaussian function similar to the
one used in serial EnSRF suggested by Hunt, pers. comm.), the optimal size of the
LEKF patches becomes larger, allowing the use of more observations. The overall LEKF
performance is then equal or better than that of the serial EnSRF with similar parameters.
We also introduced a new method to objectively estimate the optimal covariance in ation
that could be performed during the data assimilation, allowing for in ation dependent on
location and time.
In Chapter 3 the two EnKF methods as well as the three-dimensional variational
method (3DVAR) were applied to the SPEEDY primitive-equation global model (Molteni
2003). The SPEEDY model is a fast but relatively realistic model allowing a comparison of
methods that address the flrst question. Our results show that in a perfect model scenario
the EnKF greatly outperforms 3DVAR. Surprisingly, for each of the methods the 2-day
forecast ?errors of the day? are very similar to the analysis errors, but they are not similar
among difierent methods except in regions identifled as having low ensemble dimension (E-
dim), where they are also similar to bred vectors. Our results suggest that for the SPEEDY
model serial EnSRF outperforms LEKF, but that their difierence is substantially reduced if
weeitherlocalizetheLEKFerrorcovarianceorincreasetheensemblesize. Usingstochastic
perturbations further improved the performance of the LEKF. Since LEKF is much more
e?cient than serial EnSRF when using parallel computers and many observations, LEKF
(or the more e?cient algorithm LETKF, Hunt et al, 2005) would be the only feasible
choices in an operational setup.
In Chapter 4 we remove the perfect model assumption and investigate the second
question, using the NCEP/NCAR reanalysis (Kalnay et al. 1996) as the ?nature? run.
The advantage of EnKF with respect to 3DVAR is greatly reduced. When we apply
174
the model bias estimation proposed by Dee and da Silva (1998), we flnd that the full
dimensional model bias estimation fails. However, if instead we assume that the bias is
low dimensional, we obtain a substantial improvement in the EnKF analysis.
5.2 Future directions
Since the present research aims to develop a path towards an operational application of
EnKF, we would like to conclude this thesis by describing possible future directions.
Possible further investigations relevant in operations are as follows:
? The relative advantages and disadvantages between EnKF and the four-dimensional
variational method (4DVAR)
? Application of four-dimensional EnKF (4D-EnKF, Hunt et al. 2004) to treat asyn-
optic observations
? The efiect of using a nonlinear observational operator H
? The efiect of non-diagonal components of observational error covariance R
? Application to real observations
4DVAR is the mainstream in current operational data assimilation; following the ECMWF,
several major operational weather centers including the Meteo France, the Japan Meteo-
rological Agency (JMA), the UK MetO?ce (UKMO), and the Canadian Meteorological
Centre (CMC) applied 4DVAR operationally. However, there is no research thus far com-
paring EnKF with 4DVAR. There are many asynoptic observations in practice, 4DVAR
can treat them intrinsically. Alternatively, Hunt et al. (2004) proposed 4D-EnKF, as-
suming ensemble forecasts can propagate observational increments in time. When the
175
observational operator H is nonlinear, EnKF has an advantage over variational methods
in that it does not require a linearized operator or its adjoint. In addition, LEKF can
consider non-diagonal components of the observational error covariance R with less addi-
tional computation than other methods. Since more satellite observations with nonlinear
observational operator and correlated errors were available recently, it is important to
investigate these efiects.
As for model errors, research suggests that cross-correlation between forecast error
and bias error may play an important role (H. Li, pers. comm.). Since our investigation
on the model bias estimation has limitations, more investigation is necessary. If a model
bias correction scheme using both forecast and bias ensembles shows advantages, EnKF
would be much advantageous over 4DVAR. In addition, there is another issue related to
model errors that we did not discuss so far: model error covariance. Ensemble members
are expected to represent error structures, but model error covariance could afiect them
even in the absence of model biases. This may be important with real observations.
Furthermore, it is important to improve not only the analysis but also the forecast by
considering model errors. Danforth et al. (2005) have shown improved forecasts just by
keeping a constant bias as a nudging term in the forecast model using a quasi-geostrophic
global model (Marshall and Molteni 1993) and the SPEEDY model. If we estimate model
biases adaptively using a bias estimation scheme with EnKF, we may apply the same
approach as Danforth to improve forecast skill. The better forecasts may feedback in
error representation of an ensemble, improving analysis as well. However, we need to
consider the efiect of the improved forecast in the bias estimation, since the forecast
model is difierent every time.
There are further applications of ensembles representing error structures:
176
? Application in targeted or adaptive observations
? Advantage in ensemble forecasting
? Application in spatially weighted in ation parameter
If we select observing locations adaptively (targeted or adaptive observations, e.g.
Lorenz and Emanuel 1998), we can use the information of errors represented by an ensem-
ble to select the best locations to observe. It is important to observe in large error areas,
an ensemble can suggest the appropriate locations.
The ensemble generated by EnKF should be an optimal choice for ensemble forecast-
ing. Ensemble members generated by EnKF are not positive-negative pairs, Wang et al.
(2004) have shown a centered spherical simplex ensemble provides better ensemble fore-
casts than an ensemble of positive-negative pairs with a limited ensemble size. The ensem-
ble generated by EnKF represents analysis errors, Wang and Bishop (2003) showed bred
vectors weighted by the estimated analysis errors using observational locations (masked
breeding scheme, Wang and Bishop 2003) provided better ensemble forecasting.
We could also use the information of E-dimension to weight the covariance in ation
parameter in space, as suggested by Kalnay (pers. comm.). Low E-dimension areas show
dynamically important areas with good agreement of error flelds generated separately. We
can make the covariance in ation smaller in those areas, since we can trust the ensemble
more in these low E-dimensional areas.
A big challenge in operations is in applying EnKF in regional/mesoscale models.
Operationally, regional/mesoscale models provide important weather products, more di-
rectly connected to peoples lives. It is important that we apply better methods in re-
gional/mesoscale prediction systems. There are not many results on ensemble forecasting
177
or EnKF on mesoscale models, and few ensemble prediction systems are operational. Since
the regional models require lateral boundary conditions, their treatment introduces a new
problem. A global ensemble generated by EnKF may provide a good set of lateral bound-
ary conditions for regional ensemble members. In addition, the scales of weather phenom-
ena expressed by regional/mesoscale models are difierent from global models. Asynop-
tic observations capturing such mesoscale phenomena are more important and should be
treated appropriately. Meso-EnKF could be an important possible direction in operations.
178
Appendix A
The model-independent core modules of serial EnSRF and LEKF
A.1 Introduction
Since the methods of serial EnSRF and LEKF are general, they can be applied to any dy-
namical system with observations, including the simple Lorenz-96 model and the realistic
SPEEDY model. Thus, the core parts are developed in a model-independent way so that
they can be tested using a simple model and applied to a realistic model. This appendix
describes details in the model-independent modules that are written in Fortran90.
A.2 Serial EnSRF
Three subroutines are developed for serial EnSRF:
1. enkf_anal0(nx,m,hdxf,r,dxf,gain)
2. enkf_serial(nx,m,hdxf,r,dxf,gain,dxa)
3. enkf_schur(scale,dist,factor)
The flrst subroutine enkf_anal0 computes the Kalman gain using eq.(2.38), i.e.,
K = Ef(HEf)T[HEf(HEf)T +R]?1 (A.1)
As inputs, Ef, HEf and R are given, and K is returned as an output. Table A.1 gives
all the arguments of the subroutine enkf_anal0.
179
Argument IN/OUT Entity
nx IN dimension of the system
m IN ensemble size
hdfx(m) IN HEf
r IN observational error variance
dxf(nx,m) IN forecast ensemble perturbations Ef
gain(nx) OUT Kalman gain K
Table A.1: Arguments of the subroutine enkf anal0 that computes Kalman gain.
The second subroutine enkf_serial computes analysis ensemble perturbations us-
ing eqs.(2.45) and (2.48), that is,
Ea = (I?fiKH)Ef (A.2)
= Ef ?fiKHEf (A.3)
Ef, HEf, R and K are given as inputs, and Ea is returned as an output. Table A.2 gives
all the arguments of the subroutine enkf_serial.
The third subroutine enkf_schur computes the factor S(r) of the Schur product
using eq.(2.49). Table A.3 gives all the arguments of the subroutine enkf_schur.
The Schur product is applied just after computing K by enkf_anal0, and the input
of K for the subroutine enkf_serial is localized. Analysis equation eq.(2.1), that is,
?xa = ?xf +Kd (A.4)
is solved w.r.t. the ensemble mean state, where the localized Kalman gain K is applied to
avoid sampling errors in distant points. Since there is no matrix inversion or square root
computation in the above equations, no subroutine for eigenvalue decomposition or other
180
Argument IN/OUT Entity
nx IN dimension of the system
m IN ensemble size
hdfx(m) IN HEf
r IN observational error variance
dxf(nx,m) IN ensemble perturbations Ef
gain(nx) OUT Kalman gain K
dxa(nx,m) OUT analysis ensemble perturbations Ea
Table A.2: Arguments of the subroutine enkf serial that computes analysis ensemble per-
turbations.
Argument IN/OUT Entity
scale IN length scale lc
dist IN distance d
factor OUT Schur product factor S(r)
Table A.3: Arguments of the subroutine enkf schur that computes the factor of the Schur
product for localization.
181
Argument IN/OUT Entity
h(nobsl,ndiml) IN observational operator in a local patch
rdiag(nobsl) IN observational error variance in a local patch
dep(nobsl) IN observational increment d in a local patch
xb(ndiml,m) IN forecast ensemble Ef + ?xf in a local patch
xa(ndiml,m) OUT analysis ensemble perturbation Ea in a local patch
xabar(ndiml) OUT analysis ensemble mean ?xa
Table A.4: Arguments of the subroutine lekf core. ndiml, nobsl, and m denote the dimen-
sion of the local patch, the number of observations in the local patch, and ensemble size,
respectively.
matrix-related procedures is required for this EnKF module.
A.3 LEKF
A core subroutine lekf_core(h,rdiag,dep,xb,xa,xabar) is developed for LEKF. Since
the analysis equations are solved in local patches, the dimension of the local patch (ndiml)
and the number of observations in the local patch (nobsl), which are deflned as global
variables of the module, need to be determined in advance. Table A.4 gives arguments of
the subroutine lekf_core, where m denotes ensemble size.
Inside the subroutine, subroutines for real-symmetric matrix inverse, square root,
and eivenvalue decomposition are called. Each needs to be prepared separately.
182
Appendix B
Recursive fllter technique (Purser et al. 2003a)
B.1 Introduction
As introduced in Section 3.2.2, recursive fllter (RF) is applied to realize the Gaussian
shape of the horizontal correlation in 3DVAR. RF is a quasi-Gaussian fllter which has an
important advantage in computational cost.
3DVAR requires forward and adjoint transformations. Thus, given Gaussian trans-
formation G, the Cholesky decomposition
G = LLT (B.1)
provides a solution that can be used in 3DVAR. Here, L is a lower triangular matrix which
gives a forward transformation, andLT gives its adjoint. Although this method is currently
used in some operational 3DVAR systems, this matrix application is computationally
expensive.
Alternatively, RF is much cheaper. RF consists of a set of forward and backward
processes mutually adjoint:
Bi = flAi +
nX
j=1
fijBi?j (B.2)
Ci = flBi +
nX
j=1
fijCi+j (B.3)
which are nth-order RF (nth-RF), cf. eqs.(3.7) and (3.8) in the case of 4th-order. Here,
A is an input signal, and C is the flltered output signal. The subscript i denotes the space
index of grid points in one direction. Eq.(B.2) gives forward transformation, and eq.(B.3)
183
gives its adjoint. The key to implement RF is to flnd coe?cients fi and fl. This appendix
introduces RF theory based on Purser et al. (2003a) and describes details on how to flnd
the RF coe?cients especially in the case of 4th-RF.
B.2 nth-order recursive fllter
Deflne a difierential operator D(n) as
D(n) = 1?  
2?x2
2
d2
dx2 +
1
2!
?
 2?x2
2
d2
dx2
!2
+???+ 1n!
?
? 
2?x2
2
d2
dx2
!n
(B.4)
where  and ?x denote a length scale measured in numbers of grids and a grid spacing,
respectively, i.e.  ?x denotes a length scale in a physical unit. D(1) gives
D(1) = exp
?
? 
2?x2
2
d2
dx2
!
(B.5)
Let A(x) and C(x) be input and output signals of D?1(1), respectively. Consider their
spectral representations ?A(k) and ?C(k):
A(x) = 12?
Z
?A(k)eikxdk (B.6)
C(x) = 12?
Z
?C(k)eikxdk (B.7)
Then,
D(1)C(x) = 12?
Z
?C(k)exp
?
? 
2?x2
2
d2
dx2
!
eikxdk (B.8)
= 12?
Z
?C(k)exp
?
 2?x2
2 k
2
!
eikxdk (B.9)
Thus, the spectral representation of the pseudo-difierential operator D(1) becomes
?D(1) = exp
?
 2?x2
2 k
2
!
(B.10)
For an input signal with the form of delta function
A(x) = ?(x) = 12?
Z
eikxdk (B.11)
184
the transformation by ?D?1(1) is written as
?D?1(1)?(x) = 1
2?
Z
exp
?
? 
2?x2
2 k
2
!
eikxdk (B.12)
= 12? exp
?
? x
2
2 2?x2
!
(B.13)
which is nothing but a Gaussian function with the length scale  ?x. Thus, the trans-
formation ?D?1(1) gives a Gaussian fllter, and ?D?1(n) is the nth-order approximation, that
is,
C = D?1(n)A (B.14)
gives an nth-order quasi-Gaussian fllter that transforms the input signal A to get the
smoothed output C with the Gaussian shape with the length scale  ?x.
Deflne a flnite difierence operator K as
K?i = ??i?1 +2??i+1 (B.15)
where ? denotes an arbitrary variable with one-dimensional space index i (cf. eq.(3.1) in
Purser et al. 2003a). Second order derivative with flnite difierencing is written as
d2
dx2?i =
?i?1 ?2?i +?i+1
?x2 (B.16)
= ? K?x2?i (B.17)
where uniform grid spacing is assumed. Using eq.(B.17), we can rewrite eq.(B.4) as
D(n) = 1+  
2K
2 +
1
2!
?
 2K
2
!2
+???+ 1n!
?
 2K
2
!n
(B.18)
In fact, eq.(B.18) is not as accurate as what can be obtained using spectral representation,
that is given by eq.(3.12) in Purser et al. (2003a). As shown in Purser et al. (2003a), small
corrections to the coe?cients of the terms Kj(j = 1,???,n) give better results. In either
185
way, D(n) is written in the form of an nth-degree polynomial w.r.t. K. The polynomial
can be written in a difierent form as
D(n) =
nY
j=1
?
1? K?
j
!
(B.19)
cf. eq.(A.1) in Purser et al. (2003a).
Deflne a shift operator Z as
Z?i = ?i+1 (B.20)
cf. eq.(A.2) in Purser et al. (2003a). Then, eq.(B.15) can be written as
K = ?Z?1 +2?Z (B.21)
cf. eq.(A.3) in Purser et al. (2003a). Each factor of eq.(B.19) can be written as
1? K?
j
=
?
1??jZ?1
1??j
!?
1??jZ
1??j
!
(B.22)
cf. eq.(A.5) in Purser et al. (2003a), where ?j satisfles
?2j ?(2??j)?j +1 = 0 (B.23)
Thus, the quasi-Gaussian transformation (eq.(B.14)) is rewritten as
Ci =
nY
j=1
?
1??j
1??jZ?1
!?
1??j
1??jZ
!
Ai (B.24)
We can separate the transformation process (eq.(B.24)) as follows:
Bi =
nY
j=1
?
1??j
1??jZ?1
!
Ai (B.25)
Ci =
nY
j=1
?
1??j
1??jZ
!
Bi (B.26)
which are equivalent to
Bi = flAi +
nX
j=1
fijZ?jBi (B.27)
Ci = flBi +
nX
j=1
fijZjCi (B.28)
186
where
fi1 = ?1 +?2 +???+?n (B.29)
fi2 = ?1?2 +?1?3 +???+?n?1?n (B.30)
fi3 = ?1?2?3 +?1?2?4 +???+?n?2?n?1?n (B.31)
...
fin = ?1?2?3????n (B.32)
fl = (1??1)(1??2)???(1??n)
= 1?fi1 ?fi2 ?????fin (B.33)
Eqs. (B.27) and (B.28) are nothing but forward and backward processes of the nth-RF
(eqs.(refeq:nrfi2) and (refeq:nrfb2)). Thus, solving the quadratic equation (B.23) for ?j,
we get coe?cients for the nth-RF using eqs.(B.29)-(B.33). To solve the quadratic equation
(B.23), ?j are required. Since the polynomials (B.18) and (B.19) should be identical, we
get ?j by comparing the coe?cients. The polynomial (B.18) gives the numerical values
of the coe?cients according to the length scale  . The other polynomial (B.19) gives the
coe?cients
nY
j=1
?
1? K?
j
!
= 1 +
 
? 1?
1
? 1?
2
????? 1?
n
?
K
+
 1
?1?2 +???+
1
?n?1?n
?
K2
...
+ (?1)n
 1
?1?2????n
?
Kn (B.34)
In the nth-RF, n roots of an nth-degree equation give ?j. The details to get ?j are
described in the next section in the case of the 4th-RF.
187
B.3 4th-order recursive fllter
Since 4th-RF has been implemented in the 3DVAR system on the SPEEDY model, this
section describes speciflc solutions on the 4th-RF. Using the corrected polynomial given
by eq.(3.12) in Purser et al. (2003a), eq.(B.18) in the 4th-order case is written as
D(4) = 1 +  
2
2 K +
?
 4
8 +
 2
24
!
K2 +
?
 6
48 +
 4
48 +
 2
180
!
K3
+
?
 8
384 +
 6
192 +
 4
1920 +
 2
1120
!
K4 (B.35)
Alternatively, eq.(B.34) in the 4th-order case gives
D(4) =
4Y
j=1
?
1? K?
j
!
= 1 +
 
? 1?
1
? 1?
2
? 1?
3
? 1?
4
?
K
+
 1
?1?2 +
1
?1?3 +
1
?1?4 +
1
?2?3 +
1
?2?4 +
1
?3?4
?
K2
+
 
? 1?
1?2?3
? 1?
1?2?4
? 1?
1?3?4
? 1?
2?3?4
?
K3
+
 1
?1?2?3?4
?
K4 (B.36)
From the eqs.(B.35) and (B.36),
a =  
2
2 = ?
1
?1 ?
1
?2 ?
1
?3 ?
1
?4 (B.37)
b =  
4
8 +
 2
24 =
1
?1?2 +
1
?1?3 +
1
?1?4 +
1
?2?3 +
1
?2?4 +
1
?3?4 (B.38)
c =  
6
48 +
 4
48 +
 2
180 = ?
1
?1?2?3 ?
1
?1?2?4 ?
1
?1?3?4 ?
1
?2?3?4 (B.39)
d =  
8
384 +
 6
192 +
 4
1920 +
 2
1120 =
1
?1?2?3?4 (B.40)
Since eqs.(B.37)-(B.40) yield
 
x? 1?
1
? 
x? 1?
2
? 
x? 1?
3
? 
x? 1?
4
?
= x4 +ax3 +bx2 +cx+d = 0 (B.41)
the four complex roots of the 4th-degree equation (B.41) give 1/?j.
Provided ?j, ?j are computed by solving the quadratic equation (B.23):
?j = 1? ?j2 ?
s 
1? ?j2
?2
?1 (B.42)
188
As with ?j, the ones with the smaller absolute value are chosen. From eqs.(B.29)-(B.33),
coe?cients for the 4th-RF is given as
fi1 = ?1 +?2 +?3 +?4 (B.43)
fi2 = ?1?2 +?1?3 +?1?4 +?2?3 +?2?4 +?3?4 (B.44)
fi3 = ?1?2?3 +?1?2?4 +?1?3?4 +?2?3?4 (B.45)
fi4 = ?1?2?3?4 (B.46)
fl = 1?fi1 ?fi2 ?fi3 ?fi4 (B.47)
all of which are real numbers.
189
BIBLIOGRAPHY
[1] Anderson, J. L. and S. L. Anderson, 1999: A Monte Carlo Implementation of the
Nonlinear Filtering Problem to Produce Ensemble Assimilations and Forecasts. Mon.
Wea. Rev., 127, 2741-2758.
[2] Anderson, J.L., 2001: AnEnsembleAdjustmentKalmanFilterforDataAssimilation.
Mon. Wea. Rev., 129, 2884-2903.
[3] Andrews, A., 1968: A square root formulation of the Kalman covariance equations.
AIAA J., 6, 1165-1168.
[4] Atlas, R., 1997: Atmospheric observations and experiments to assess their usefulness
in data assimilation. J. Meteor. Soc. Japan, 75(1B), 111-130.
[5] Barker, D. M., W. Huang, Y.-R. Guo, A. J. Bourgeois, and Q. N. Xiao, 2004: A
Three-Dimensional Variational Data Assimilation System for MM5: Implementation
and Initial Results. Mon. Wea. Rev., 132, 897-914.
[6] Bouttier, F., 1996: The Kalman Filter. Proc. Seminar on Predictability, Vol.1.,
ECMWF, Reading, Berkshire, UK, 221-245.
[7] Bouttier, F. and P. Courtier, 1999: Data assimilation concepts and methods. Meteo-
rological Training Course Lecture Series, ECMWF, 75pp.
[8] Burgers, G., P. J. van Leeuwen, and G. Evensen, 1998: Analysis Scheme in the
Ensemble Kalman Filter. Mon. Wea. Rev., 126, 1719-1724.
[9] Bishop, C. H., B. J. Etherton, and S. J. Majumdar, 2001: Adaptive Sampling with
Ensemble Transform Kalman Filter. Part I: Theoretical Aspects. Mon. Wea. Rev.,
129, 420-436.
190
[10] Cane, M. A., A. Kaplan, R. N. Miller, B. Tang, E. C. Hackert, and A. J. Busalacchi,
1996: Mapping tropical Paciflc sea level: Data assimilation via a reduced state space
Kalman fllter. J. Geophys. Res., 101(C10), 22,599-22,617.
[11] Corazza, M., E. Kalnay, D. J. Patil, E. Ott, J. Yorke, I. Szunyogh, M. Cai, 2002: Use
of the breeding technique in the estimation of the background error covariance matrix
for a quasigeostrophic model. AMS Symposium on Observations, Data Assimilation
and Probabilistic Prediction, Orland, Florida, 154-157.
[12] Corazza, M., E. Kalnay, D. J. Patil, I. Szunyogh, and S.-C. Yang, 2003: Keeping the
bred vectors young: Impact on data assimilation. The presentation flle available at
http://www.atmos.umd.edu/ ekalnay/.
[13] Daley, R., 1991: Atmospheric Data Analysis. Cambridge University Press, 457pp.
[14] Danforth, C., E. Kalnay, and T. Miyoshi, 2005: Estimating and Correcting Global
Weather Model Error. In preparation.
[15] Dee, D. P., 1995: On-line Estimation of Error Covariance Parameters for Atmospheric
Data Assimilation. Mon. Wea. Rev., 123, 1128-1145.
[16] Dee, D. P. and A. M. da Silva, 1998: Data assimilation in the presence of forecast
bias. Quart. J. Roy. Meteor. Soc., 126, 269-295.
[17] Ehrendorfer, M. and R. M. Errico, 1995: Mesoscale Predictability and the Spectrum
of Optimal Perturbations. J. Atmos. Sci., 52, 3475-3500.
[18] Etherton, B. J. and C. H. Bishop, 2004: Resilience of Hybrid Ensemble/3DVAR
Analysis Schemes to Model Error and Ensemble Covariance Error. Mon. Wea. Rev.,
132, 1065-1080.
191
[19] Evensen, G., 1994: Sequential data assimilation with a nonlinear quasi-geostrophic
model using Monte Carlo methods to forecast error statistics. J. Geophys. Res., 99
(C5), 10143-10162.
[20] Fukumori I. and P. Malanotte-Rizzoli, 1995: An approximate Kalman fllter for ocean
data assimilation: An example with an idealized Guls Stream model. J. Geophys.
Res., 100 (C4), 6777-6793.
[21] Gaspari, G. and S. E. Cohn, 1999: Construction of correlation functions in two and
three dimensions. Quart. J. Roy. Meteor. Soc., 125, 723-757.
[22] Gelb, A., J. F. Kasper, R. A. Nash, C. F. Price, and A. A. Sutherland, 1974: Applied
Optimal Estimation. The M.I.T. Press, 374pp.
[23] Hamill, T. M. and C. Snyder, 2000: A Hybrid Ensemble Kalman Filter-3D Variational
Analysis Scheme. Mon. Wea. Rev., 128, 2905-2919.
[24] Hamill, T. M., J. S. Whitaker, and C. Snyder, 2001: Distance-Dependent Filtering of
Background Error Covariance Estimates in an Ensemble Kalman Filter. Mon. Wea.
Rev., 129, 2776-2790.
[25] Hamill, T. M., 2003: Ensemble-based Atmospheric Data Assimilation: A Tutorial.
University of Colorado and NOAA-CIRES Climate Diagnostics Center, Boulder, Col-
orado, 46pp.
[26] Houtekamer, P. L. and H. L. Mitchell, 1998: Data Assimilation Using an Ensemble
Kalman Filter Technique. Mon. Wea. Rev., 126, 796-811.
[27] Houtekamer, P. L. and H. L. Mitchell, 1999: Reply. Mon. Wea. Rev., 127, 1378-1379.
192
[28] Houtekamer, P. L. and H. L. Mitchell, 2001: A Sequential Ensemble Kalman Filter
for Atmospheric Data Assimilation. Mon. Wea. Rev., 129, 123-137.
[29] Houtekamer, P. L. and H. L. Mitchell, 2003: Practical Ensemble Data Assimilation.
Seminar on Recent developments in data assimilation for atmosphere and ocean, 8 to
12 September 2003, ECMWF, 193-210.
[30] Houtekamer, P. L., H. L. Mitchell, G. Pellerin, M. Buehner, M. Charron, L. Spacek,
and B. Hansen, 2005: Atmospheric Data Assimilation with an Ensemble Kalman
Filter: Results with Real Observations. Mon. Wea. Rev., 133, 604-620.
[31] Hunt, B. R., E. Kalnay, E. J. Kostelich, E. Ott, D. J. Patil, T. Sauer., I. Szunyogh, J.
A. Yorke and A. V. Zimin, 2004: Four-dimensional ensemble Kalman flltering. Tellus,
56A, 273-277.
[32] Ide, K., P. Courtier, M. Ghil, and A. C. Lorenc, 1997: Unifled notation for data
assimilation: Operational, sequential, and variational. J. Meteor. Soc. Japan, 75
(1B), 181-189.
[33] Jazwinski, A. H., 1970: Stochastic Processes and Filtering Theory. Academic Press,
376pp.
[34] Kalman, R. E., 1960: A New Approach to Linear Filtering and Prediction Problems.
J. Basic Eng., Trans. ASME, 35-45.
[35] Kalnay, E. and Z. Toth, 1994: Removing growing errors in the analysis cycle.
Preprints, Tenth Conference on Numerical Weather Prediction, Amer. Meteor. Soc.,
212-216.
193
[36] Kalnay, E., M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell,
S. Saha, G. White, J. Woollen, Y. Zhu, A. Leetmaa, B. Reynolds, M. Chelliah, W.
Ebisuzaki, W. Higgins, J. Janowiak, K. C. Mo, C. Ropelewski, J. Wang, R. Jenne
and D. Joseph: The NCEP/NCAR 40-Year Reanalysis Project. Bull. Amer. Meteor.
Soc., 77(3), 437-471.
[37] Kalnay, E., 2003: Atmospheric modeling, data assimilation and predictability. Cam-
bridge University Press, 341pp.
[38] Kistler, R., E. Kalnay, W. Collins, S. Saha, G. White, J. Woollen, M. Chelliah, W.
Ebisuzaki, M. Kanamitsu, 2001: The NCEP-NCAR 50-Year Reanalysis: Monthly
Means CD-ROM and Documentation. Bull. Amer. Meteor. Soc., 82(2), 247-267.
[39] Lorenz, E., 1996: Predictability - A problem partly solved. Proc. Seminar on Pre-
dictability, Reading, United Kingdom, ECMWF, 1-18.
[40] Lorenz, E. and K. Emanuel, 1998: Optimal sites for supplementary weather observa-
tions. J. Atmos. Sci., 55, 399-414.
[41] Marshall, J. and F. Molteni, 1993: Toward a Dynamical Understanding of Planetary-
Scale Flow Regimes. J. Atmos. Sci., 50, 1792-1818.
[42] Miyoshi, T., 2004: Classical methods tour of advanced data assimilation using Lorenz
96 model. A flnal report submitted to METO658E, Department of Meteorology, Uni-
versity of Maryland, College Park.
[43] Miyoshi, T., 2004: A review of the deterministic ensemble Kalman flltering meth-
ods and related techniques. Masters Scholarly Paper. Department of Meteorology,
University of Maryland, College Park, 20pp.
194
[44] Miyoshi, T. and E. Kalnay, 2005: Ensemble forecasting using orthogonalized and
stochastically seeded bred vectors in the Lorenz 96 model with orography. In prepa-
ration.
[45] Molteni, F., 2003: Atmospheric simulations using a GCM with simplifled physical
parameterizations. I: model climatology and variability in multi-decadal experiments.
Clim. Dyn., 20, 175-191.
[46] Nicolis, C., 2003: Dynamics of Model Error: Some Generic Features. J. Atmos. Sci.,
60, 2208-2218.
[47] Ott, E., B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E.
Kalnay, D. J. Patil, and J. A. Yorke, 2002: Exploiting local low dimensionality of
the atmospheric dynamics for e?cient Kalman flltering. ArXiv:arc-ive/paper 020358,
http://arxiv.org/abs/physics/020358.
[48] Ott, E., B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E. Kalnay,
D. J. Patil, and J. A. Yorke, 2004: A local ensemble Kalman fllter for atmospheric
data assimilation. Tellus, 56A, 415-428.
[49] Patil, D. J., B. R. Hunt, E. Kalnay, J. A. Yorke, and E. Ott, 2001: Local Low
Dimensionality of Atmospheric Dynamics. Phys. Rev. Lett., 5878-5881.
[50] Parrish, D. F. and J. C. Derber, 1992: The National Meteorological Center?s Spectral
Statistical-Interpolation Analysis System. Mon. Wea. Rev., 120, 1747-1763.
[51] Purser, R. J., W. S. Wu, D. F. Parrish, and N. M. Roberts, 2003a: Numerical Aspects
of the Application of Recursive Filters to Variational Statistical Analysis. Part I:
195
Spatially Homogeneous and Isotropic Gaussian Covariances. Mon. Wea. Rev., 131,
1524-1535.
[52] Purser, R. J., W. S. Wu, D. F. Parrish, and N. M. Roberts, 2003b: Numerical
Aspects of the Application of Recursive Filters to Variational Statistical Analysis.
Part II: Spatially Inhomogeneous and Anisotropic General Covariances. Mon. Wea.
Rev., 131, 1536-1548.
[53] Szunyogh, I., E. J. Kostelich, G. Gyarmati, D. J. Patil, B. R. Hunt, E. Kalnay, E.
Ott, and J. A. Yorke, 2004: Assessing a local ensemble Kalman fllter: Perfect model
experiments with the NCEP global model. Tellus, in print.
[54] Tippett, M. K., J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker,
2003: Ensemble Square Root Filters. Mon. Wea. Rev., 131, 1485-1490.
[55] van Leeuwen, P. J., 1999: Comment on ?Data Assimilation Using an Ensemble
Kalman Filter Technique?. Mon. Wea. Rev., 127, 1374-1377.
[56] Whitaker, J. S. and T. M. Hamill, 2002: Ensemble Data Assimilation without Per-
turbed Observations. Mon. Wea. Rev., 130, 1913-1924.
[57] Wang, X. and C. H. Bishop, 2003: A Comparison of Breeding and Ensemble Trans-
form Kalman Filter Ensemble Forecast Schemes. J. Atmos. Sci., 60, 1140-1158.
[58] Wang, X., C. H. Bishop, S. J. Julier, 2004: Which Is Better, an Ensemble of Positive-
Negative Pairs or a Centered Spherical Simplex Ensemble?. Mon. Wea. Rev., 132,
1590-1605.
[59] Wilks, D. S., 1995: Statistical Methods in the Atmospheric Sciences, Academic Press,
467pp.
196
[60] Zupanski, D., 1997: A General Weak Constraint Applicable to Operational 4DVAR
Data Assimilation Systems. Mon. Wea. Rev., 125, 2274-2292.
197