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ﬁltering (EnKF) system. Several approaches to EnKF for atmospheric systems have been
proposed but not systematically compared. The sensitivity of EnKF to the imperfectionsof forecast models is unclear. This research explores two questions: 1. What are therelative advantages and disadvantages of two promising EnKF methods? 2. How large
are the eﬀects 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 ﬁlter (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 covarianceinﬂation.
In Chapter 3 we apply the two EnKF methods and the three-dimensional variational 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 analysis errors, but they are not similar among diﬀerent methods except in low ensemble
dimensional regions. Overall, serial EnSRF outperforms LEKF, but their diﬀerence issubstantially 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 weapply the model bias estimation proposed by Dee and da Silva (1998), we ﬁnd that thefull 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 MODELby
Takemasa Miyoshi

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulﬁllmentof the requirements for the degree ofDoctor 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

ACKNOWLEDGMENTSFirst of all, I express my deep gratitude to Prof. Eugenia Kalnay for her advisingand supports during the entire processes of my graduate study. She always supported meboth academically and mentally, which was indispensable for me to complete this hardtask.
Secondly, I thank the Japanese Government, especially the National Personnel Authority, 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 theiruseful 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 discussions 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 ﬁlter theory. Kohei Aranami kindly
volunteered to let me use his powerful personal computer remotely. Some of the analysis/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 tocomplete 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 Figuresxii

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.3

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 ﬁltering (serial EnSRF) . . . . . . . . . 18

2.2.5

Local ensemble Kalman ﬁltering (LEKF) . . . . . . . . . . . . . . . 22

2.2.6

Observational error covariance localization . . . . . . . . . . . . . . . 23

2.2.7

Covariance inﬂation . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.2.8

Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . 26

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

3

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

Data assimilation on the SPEEDY primitive-equation model
3.1

3.2

3.3

3.4

3.5

46

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.1.2

The SPEEDY model . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

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

EnKF implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.1

Serial EnSRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.3.2

LEKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.1

Observational network . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.4.2

Data assimilation experiments . . . . . . . . . . . . . . . . . . . . . 67

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5.1

3DVAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

3.5.2

Serial EnSRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
v

3.6

3.7

3.8

3.5.3

LEKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

3.5.4

Timing results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

3.5.5

Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Sensitivity experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.6.1

Response with a diﬀerent observational network

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

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 ﬁelds 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

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

4 Data assimilation in the presence of model errors
4.1

4.2

. . . . . . . . . . . 101

138

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.1.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

4.1.2

OSSE methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

4.1.3

NCEP/NCAR reanalysis (NNR) . . . . . . . . . . . . . . . . . . . . 139

Model error statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
vi

4.2.1Model error biases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

4.2.2

EOFs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

4.3

Eﬀects of model errors on data assimilation . . . . . . . . . . . . . . . . . . 146

4.4

Model bias estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

4.5

4.4.1

Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

4.4.2

Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 158

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 ﬁlter technique (Purser et al. 2003a)

183

B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
B.2 nth-order recursive ﬁlter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
B.3 4th-order recursive ﬁlter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188

Bibliography

190

vii

LIST OF TABLES

2.1

The 1-year mean analysis RMSE of each method with optimal parametersettings 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 ﬁxed 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 isapplied, and the local patch parameter is ﬁxed 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ﬂation parameter is applied. l2 is ﬁxed to l2 = 0. In the cases with observational error covariance localization, l is ﬁxed 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 ﬁeld in meters using 3DVAR on theSPEEDY 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. ”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 zonallyand horizontally uniform background error standard deviation is considered, respectively. ”NO ASSIM” denotes the case that no data assimilation
is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.4

Analysis RMSE of 500hPa height ﬁeld 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 3DVARin 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 ﬁlter divergence. 4% multiplicative
covariance inﬂation is applied. . . . . . . . . . . . . . . . . . . . . . . . . . . 78

3.5

Analysis RMSE of 500hPa height ﬁeld in meters using LEKF 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-monthspin-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, andthe local patch parameter l is chosen as integers from 1 to 6. ”D” denotes
ﬁlter divergence. 4% multiplicative covariance inﬂation is applied.

ix. . . . . 86

3.6

Analysis RMSE of 500hPa height ﬁeld in meters using LEKF on the SPEEDY
model when the sparse observational network is applied and the parameterl2 is changed. The ensemble size (NBV) is chosen to be 20 or 30, the local patch parameter is ﬁxed 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 memberson 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 theSPEEDY model when the sparse observational network is applied. The local
patch parameter l is chosen as an optimal setting shown in the parenthesis,
and l2 is ﬁxed 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 ﬁxed to 4%, the parameter l2 in LEKF is ﬁxed to 0. . . . . . . 97

3.10 Comparison of analysis RMSE of 500hPa height ﬁeld in meters on theSPEEDY model when the sparse observational network is applied. For
each ensemble size (NBV), best results are chosen for serial EnSRF andLEKF 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 assimilatedusing a parallel computer with 100 computational nodes. Parallelizationeﬃciency is assumed as 80% and 20% for LEKF and serial EnSRF, respectively. We assume the increased observations require 50 and 100 times morecomputation for LEKF and serial EnSRF, respectively.

. . . . . . . . . . . 100

3.13 Analysis RMSE of 500hPa height ﬁeld 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.

4.1. . . . . . . . . 112

Analysis RMSE of 500hPa height ﬁeld using EnKF with 10 ensemble members and low order bias estimation using various numbers of bias basic
ﬁelds. ”Uncorrected” denotes the case without bias correction. ”Meanbias”, ”EOF1”, and ”EOF2” denotes the cases with low-order bias correction using the mean bias ﬁeld, using the mean bias ﬁeld and ﬁrst EOF, andusing the mean bias ﬁeld and ﬁrst 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 thedimension 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

Deﬁnition of the local regions around a center point (C) for LEKF in theLorenz-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 thatis 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 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 σ = 5and δ = 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 ﬁxed to 1. The minimum error 0.20 is observed whenl = 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 ﬁxed 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.

. . 402.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 ﬁxed to l = 15 and l2 = 0. The horizontal axis shows thelength 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 statisticalstrength 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 forps and u are Pa and m/s, respectively. Strong spatial dependence especiallyin zonal structures is observed, but their noisiness indicates sampling errors. 55
xiii

3.3

Background error horizontal correlations. The top panel shows verticalstructure 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 thevertical axis shows latitudes. The horizontal axis shows horizontal lengthscale in grid spacing. The shades show correlation. . . . . . . . . . . . . . . 56

3.4

Vertical background error correlations of u-wind (solid lines with whitecircle), v-wind (dashed lines with black circle), temperature (short dashedlines with white square), and moisture (long-short dashed lines with blacksquare). The correlations are averaged horizontally. . . . . . . . . . . . . . . 58

3.5

Spatial response of the 3DVAR in the case of a single observation with theobservational increment of 1.0m/s at mid-latitude (top panel in m/s) andhigh-latitude (bottom panel in m/s) at the 4th level (σ = 0.51). . . . . . . . 60

3.6Inter-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 diﬀerence between truth
and ﬁrst guess (bottom panel). The ﬁelds are u-wind ﬁelds (in m/s) at the
4th level.

3.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Deﬁnition of the local regions around a center point for LEKF in theSPEEDY 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, thearea 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 percents 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 northern hemisphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.12 Schematic of data assimilation experiments on the SPEEDY model. ”DA”
denotes the data assimilation system. The ”initial” ﬁeld is the output by theSPEEDY model which is diﬀerent from the ”analysis” ﬁeld because the highfrequency components are smoothed out by the spectral transformation of
the model. The diﬀerence between ”initial” and ”truth” is deﬁned as theanalysis error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.13 Analysis RMSE of 500hPa height ﬁeld in the case of 3DVAR with thedense observational network (solid line), sparse observational network using
the same background error statistics based on the dense network (dashedline), 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). . 71

xv

3.14 Analysis RMSE of 500hPa height ﬁeld 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). . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.15 Horizontal structure of analysis RMSE of height ﬁeld in the case of 3DVARwith horizontally uniform background error standard deviation and thesparse observational network. Shades and contours show analysis RMSE
(in meters) and mean analysis ﬁelds (in meters), respectively. Dots indicate observational locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.16 Horizontal structure of analysis error bias (top panel) and standard deviation (bottom panel) of height ﬁeld (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 ﬁelds. Dots indicate observational
locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.17 Zonal structure of analysis RMSE of u-wind ﬁeld (top panel, in m/s) andtemperature ﬁeld (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 ﬁelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

xvi

3.18 Analysis RMSE of 500hPa height ﬁeld 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 Schurproduct is chosen as optimal. For comparison, the RMSE of 3DVAR isshown as the short dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.19 Horizontal structure of analysis RMSE of height ﬁeld, 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 ﬁelds (in meters). Dots indicate observational locations. . . . . . . 80
3.20 Horizontal structure of analysis error bias (top panel) and standard deviation (bottom panel) of height ﬁeld (in meters), similar to Fig.3.16, in thecase 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 ﬁelds. Dots indicate observational 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 thesparse 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 speciﬁc
humidity (q) ﬁelds, 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, theRMSE of 3DVAR is shown as the short dashed line. . . . . . . . . . . . . . 83

xvii

3.22 Zonal structure of analysis RMSE of u-wind ﬁeld (top panel, in m/s) andtemperature ﬁeld (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 ﬁelds, respectively. . . . . . . . . . . . . 85
3.23 Analysis RMSE of 500hPa height ﬁeld 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 ﬁeld, 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 ﬁelds (in meters). Dots indicate observational locations. . . . 89
3.25 Horizontal structure of analysis error bias (top panel) and standard deviation (bottom panel) of height ﬁeld (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 ﬁelds. 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 speciﬁc humidity (q)ﬁelds, respectively. The observational error standard deviations are shownas 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 ﬁeld (top panel, in m/s) andtemperature ﬁeld (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 ﬁelds, respectively. . . . . . . . . . . . . 93
3.28 Analysis RMSE of 500hPa height ﬁeld in the case of the sparse observational
network using 3DVAR (short dashed line), serial EnSRF (solid line) andLEKF (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 ﬁeld in the case of the realistic observational network using 3DVAR (short dashed line) and serial EnSRF (solidline). 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 ﬁelds (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 is30, and the localization scale is 3.0. Dots indicate the observational locations.104xix

3.31 Zonal structure of analysis RMSE of u-wind ﬁeld (top panel) and temperature ﬁeld (bottom panel) in the case of 3DVAR with the realistic observational network. Shades and contours show analysis RMSE and mean
analysis ﬁelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.32 The same ﬁgure 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 ﬁeld (a), hight ﬁeld (b), temperature ﬁeld(c) and speciﬁc humidity ﬁeld (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
3.34 Impact of the vertical correlation in the analysis RMSE of 500hPa heightﬁeld using LEKF with 10 (black lines) and 20 (blue lines) ensemble members. The black solid line and black short-dashed line show the cases withand without the vertical correlation, respectively. The blue dotted line andblue long-dashed line show the cases with and without the vertical correlation, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
3.35 The same ﬁgure as Fig.3.34 but for the full vertical correlation. . . . . . . . 111
3.36 Analysis RMSE of 500hPa height ﬁeld using LEKF with 10 ensemble members. 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. . . . . . . . . . . . . . . 113

xx

3.37 Analysis RMSE of 500hPa height ﬁeld using LEKF with 10 ensemble members. Here, observational error covariance localization is applied. Dashedand solid lines show the cases with and without random perturbation addition, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
3.38 Background error ﬁeld (shaded) and analysis increment ﬁeld (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 ﬁelds in the cases of 3DVAR, serial
EnSRF with 30 ensemble members and the localization scale of 3, andLEKF with 30 ensemble members, the local patch parameter l = 3 andcovariance localization scale σ = 1.5. . . . . . . . . . . . . . . . . . . . . . . 119
3.40 The growth of the initial analysis errors and ensemble perturbations in thecases 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 ﬁelds (ananlysis/forecast minus truth; shaded) and the ensemble
spread of 500hPa height ﬁeld (contour) at an arbitrary time. The top andbottom panels show analysis and background (6-hour forecast), respectively,
both at the same time. Some parts have similar structures indicated byrectangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

xxi

3.42 2-day forecast error ﬁeld (shaded) and analysis error ﬁeld (contour) of uwind at the 4th vertical level (σ = 0.51) at an arbitrary time in the caseof serial EnSRF with 10 ensemble members and the localization scale of 1.
The correlation between the two ﬁelds 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
ﬁelds are 0.79. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
3.44 The same as Fig.3.42, but in the case of 3DVAR. The correlation between
the two ﬁelds are 0.74. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
3.45 Temporal mean of correlation coeﬃcients between 2-day forecast error andanalysis error ﬁelds 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-windat the 4th level all valid at the same time chosen arbitrarily, in the cases ofserial 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 anarbitrary time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
3.52 Scatter plot of E-dimension and explained variance at the 4th level at anarbitrary 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 ﬁelds.

4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

SPEEDY model biases of q, z [m] and ps [Pa] relative to the NCEP/NCAR
reanalysis ﬁelds. Each vertical level has diﬀerent scales for q ﬁelds, it isomitted. 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 ﬁelds at the bottom level(shades) and surface pressure ﬁelds (contour) of the ﬁrst and second EOFs,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

4.4Time series of expansion coeﬃcients of the ﬁrst to fourth EOFs for a week.
Solid, dashed, dotted, dash-dotted lines show ﬁrst, second, third, and fourth
EOFs, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

4.5Time series of expansion coeﬃcients of the ﬁrst EOF (top panel) and thethird (solid line) and fourth (dotted line) EOFs (bottom panel) for 50 days. 149

xxiii

4.6

Analysis RMSE of 500hPa height ﬁeld 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 ﬁeld of 3DVAR and EnKF with 30 ensemble members in the case of the sparse and realistic observational networks inthe 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. . . . . . . . . . . . . . . . . . . . . . . . . 151

4.8

Analysis RMSE at the all pressure levels temporally averaged for one monthafter 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 realistic observational network, respectively. The four panels (a), (b), (c) and(d) correspond to u-wind, height, temperature (T ) and speciﬁc humidity (q)ﬁelds, respectively. The observational error standard deviations are shownas thin solid lines wherever applicable. . . . . . . . . . . . . . . . . . . . . . 153

4.9

Analysis RMSE of 500hPa height ﬁeld of EnKF with 10 ensemble membersin the case of the sparse observational networks in the presence of model
errors. The dashed and solid lines show the cases with and without Deeand da Silva’s full dimensional bias estimation, respectively. . . . . . . . . . 159

xxiv

4.10 Analysis RMSE of 500hPa height ﬁeld for EnKF with low order bias estimation (green short-dashed line), EnKF with constant bias subtraction(black dotted line), 3DVAR with constant bias subtraction (blue dottedline). Solid lines show the cases without bias estimation in the cases of
EnKF (black line) and 3DVAR (blue line), the same as solid and longdashed lines of Fig.4.6. Thin red line shows 6-hour forecast errors initiatedfrom the reanalysis ﬁelds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.11 Values of all the seven components of the bias variable when the mean bias
ﬁeld 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 ﬁelds (shades inmeters) in the case of imperfect model experiments using EnKF with 10
ensemble members. Contours show mean analysis ﬁelds (in meters). Top
and bottom panels show the cases without and with bias estimation, respectively. Dots indicate the observational locations. . . . . . . . . . . . . . 164
4.13 Analysis error bias ﬁelds of 500hPa height (shades in meters) in the caseof imperfect model experiments using EnKF with 10 ensemble members.
Contours show mean analysis ﬁelds (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 ﬁelds of 500hPa height (in meters) inthe case of imperfect model experiments using EnKF with 10 ensemble
members. Contours show mean analysis ﬁelds (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 monthafter the initial one-month spin-up period (112 samples) in the cases of
EnKF with low order bias estimation (green short-dashed line), EnKF withconstant 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), thesame 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 speciﬁc humidity (q) ﬁelds, 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 ﬁeld (top panel in m/s) andtemperature ﬁeld (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
ﬁelds, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
4.17 The same ﬁgure as Fig.4.16 but with bias estimation. . . . . . . . . . . . . . 170

xxvi

Chapter 1
Introduction
1.1 BackgroundIn 2004, ten typhoons hit Japan, bringing tremendous damages to the local economy andpeople’s lives. This was the worst typhoon season in 50 years. Better preparation for suchextreme 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 thatare 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 initial 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 handto 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 constantpressure lines, constant temperature lines, etc. However, even among well-trained meteorologists, the analyzed charts are diﬀerent. For high quality NWP, it is essential togenerate an ”objective analysis” with less uncertainty given the same amount of information.
Since the initial conditions are crucial to NWP, the objective analysis (data assimilation) 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 the1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 aresult, the eﬀect of errors in the initial condition is more important. The relative impact of introducing a new data assimilation system has increased, and data assimilationis 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 variational method (4DVAR) and Kalman ﬁltering (KF) have been explored and some havealready 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 noobservational instrument without errors. Even if we have a perfect instrument, a grid pointvalue for an NWP model variable will be diﬀerent from the perfect observation because
the grid point value represents a mean state for the surrounding area. The diﬀerence 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, nomodel that simulates Mother Nature can be perfect. In addition, the atmospheric system
is nonlinear, which introduces dynamically chaotic behavior. Even if we have a perfectmodel, 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 were2

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 eﬀects 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 thelarge number of degrees of freedom of NWP models. Even if we know a best methodfor 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 bemade in atmospheric data assimilation. The simpliﬁcations will decrease as computational
capability increases.
Kalman ﬁltering (KF, Kalman 1960) is a theory of data assimilation that optimallycombines 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 simpliﬁcations.
In the ensemble formulation of Kalman ﬁltering, a.k.a. ensemble Kalman ﬁltering (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 ﬁlter (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 to3

outperform the PO method.
Several SRF methods have been proposed: an ensemble transform Kalman ﬁlter(ETKF, Bishop et al. 2001), an ensemble adjustment Kalman ﬁlter (EAKF, Anderson
2001), a serial EnSRF (Whitaker and Hamill 2002), all of which are eﬀective when observational 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 diﬀerent methodsof SRF. In addition, as Tippett et al. (2003) mentioned, the formulation of SRFs is nottheoretically unique, and it is not clear if there is any optimal choice among diﬀerent
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 ashort time interval for the data assimilation cycle if the observations are frequent enough sothat 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 totreat model biases in data assimilation. Assuming that model biases are another set ofprognostic 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 forecastanalysis system. However, as mentioned in the previous section, at least four ways toimplement EnKF on NWP models have been suggested, but the relative advantages anddisadvantages among them are not clear. In addition, EnKF has been tested most extensively 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 twoquestions:
1. What are the relative advantages of diﬀerent implementations of EnKF?
2. What is the eﬀect 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 canbe divided into two groups: serial treatment of observations and simultaneous treatmentof observations. Thus, two most eﬃcient methods of each group, a serial EnSRF (Whitakerand Hamill 2002) and LEKF (Ott et al. 2002; 2004), are chosen. The present researchaims to investigate the following questions:
1. What are the relative advantages and disadvantages between a serial EnSRF andthe LEKF?
2. How large are the eﬀects 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 experiments on a primitive-equation global model are performed. The model used in thepresent research is the SPEEDY model (Simpliﬁed Parameterizations, primitivE-Equation
DYnamics, by Molteni 2003). The SPEEDY model has simpliﬁed physical processes, theresolution 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 ), speciﬁc humidity (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 toobservational errors. To investigate the ﬁrst question, we ﬁrst 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 representativeof 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 theparameters are adjusted. Three types of observational locations are tested; a regularly
distributed dense network, a regularly distributed sparse network, and an irregularly distributed 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 et6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 modeland 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
ﬁelds 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 eﬀects 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 Lorenz96 model (Lorenz 1996; Lorenz and Emanuel 1998). Section 2.2 introduces the theoryof 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 diﬀerence 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 perfectmodel assumption using the simple Lorenz-96 model. This gives an idea of how eachmethod works and how it is diﬀerent from each other. In addition, the results ensure theprogram developed for the present research is correctly coded by comparing the performance with those of Whitaker and Hamill (2002) and Ott et al. (2004). In this section, the7three methods (3DVAR/OI, serial EnKF and LEKF) and full KF are implemented andtested on the Lorenz-96 model, the methods are compared to each other. Here, full KF
means that all ﬁve KF equations, introduced in Section 2.2, are solved explicitly without
the low-rank approximation that is essential in EnKF.
Chapter 3 discusses the ﬁrst question under the perfect model assumption. In thischapter, 3DVAR, serial EnKF and LEKF are introduced separately, and implementationsof each method are described. After descriptions of implementation, the results are shownand compared. Sensitivities to diﬀerent parameters, observational networks, and localization types are discussed. The characteristics of the analysis and forecast errors arealso discussed. For latter experiments, one EnKF method (serial EnSRF) is chosen inthis 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 reanalysis ﬁelds as the nature run, we introduce model errors. First, the model error statistics
are computed, accumulated by the diﬀerence between the reanalysis ﬁelds and the model
forecasts initiated by the reanalysis ﬁelds. Then, data assimilation experiments assuming
the reanalysis ﬁelds as truth are performed without model bias estimation, which clariﬁes
the eﬀects 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 2Theory 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 diﬀerent from each other. Here, we develop the core modules used for thisentire research. The results shown in this chapter support the correctness of the Fortranprograms.

2.2 Theoretical review
This section describes a review on theory of data assimilation methods that are developed in the present research. Some basics in estimation theory are brieﬂy introduced sothat the present dissertation is self-contained. For more comprehensive introduction anddetails 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 ﬁltering, Jazwinski (1970) provides amore mathematically precise introduction and discussion. Gelb et al. (1974) also provide
a more comprehensive and practical introduction on Kalman ﬁltering. Bouttier (1996)
introduces Kalman ﬁltering 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 ananalysis 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 diﬀerent
spaces. H converts prognostic variables x to the space of observational data yo and iscalled an ”observational operator”. H is not necessarily a linear mapping.
Deﬁne diﬀerences from the true state xt :
δxa = xa − xt

(2.2)

= xf − xt(2.3)

δyo = yo − H(xt )

(2.4)

δxf

as residual errors of analysis, forecast, and observation, respectively. Subtracting xt fromthe 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 isaccurate as long as δxf is small. The analysis error covariance Pa is written as
D

Pa = δxa (δxa )>

E

(2.6)

where h·i denotes the statistical expected value. Substituting eq.(2.5) into eq.(2.6), we
10

note
a

P

¿³

=

δx + K(δy − Hδx )
f

o

f

´³

D

=

´> À

δx + K(δy − Hδx )
f

o

E

D

f

(I − KH)δxf (δxf )> (I − KH)> + 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 deﬁned
as
D

Pf

δxf (δxf )>

=
D

R =

δyo (δyo )>

E

(2.8)

E

(2.9)

and the forecast error and observational error are assumed to be uncorrelated:DE

δxf (δyo )> = 0(2.10)

We assume the error probability distribution function (PDF) is Gaussian. Then the maximum likelihood state minimizes the total analysis error variance. That is, the optimal
weight gives the minimum of the trace of Pa . Thus,∂(trace(Pa )) = 0∂K(2.11)

Using the formulas´∂ ³
trace(XYX> ) = X(Y + Y> )∂X∂
(trace(XY)) = Y>∂X

(2.12)(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)Pf H> + 2KR = 0
11(2.14)

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 = Pf H> (HPf H> + R)−1

(2.15)For simplicity, deﬁne S as
S = HPf H> + R(2.16)

Substituting eq.(2.15) into eq.(2.7), we get the analysis error covariance
Pa = (I − Pf H> S−1 H)Pf (I − Pf H> S−1 H)> + Pf H> S−1 R(Pf H> S−1 )>= Pf − 2Pf H> S−1 HPf + Pf H> S−1 HPf H> S−1 HPf + Pf H> S−1 RS−1 HPf= Pf − 2Pf H> S−1 HPf + Pf H> S−1 (HPf H> + R)S−1 HPf= Pf − 2Pf H> S−1 HPf + Pf H> S−1 SS−1 HPf= Pf − Pf H> (HPf H> + R)−1 HPf= (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)

Deﬁne the mapping that gives the true time evolution to be M t . Then
xti+1 = M t (xti )12(2.19)

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 = 0Dηη >(2.21)

E= Q(2.22)

Then, the forecast error covariance matrix is written as
Pfi

Dδxfi (δxfi )>

=

E

D

(xfi − xti )(xfi − xti )>

=

¿³

=
¿³

=
'

D¡

E
´³

M (xai−1 ) − M t (xti−1 )M (xai−1 ) − M (xti−1 ) + ηMδxai−1 + η

¢¡

´> À

M (xai−1 ) − M t (xti−1 )´³Mδxai−1 + η

M (xai−1 ) − M (xti−1 ) + η

¢> E

D

Mδxai−1 (Mδxai−1 )> + Mδxai−1 η > + η(Mδxai−1 )> + ηη >

=
D

=

´> À

E

D

Mδxai−1 (δxai−1 )> M> + ηη >E

E= MPai−1 M> + 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 isvalid as long as δxa is small. In the derivation of eq.(2.23), the analysis error δxa and therandom error η are assumed to be uncorrelated, that is,DEδxa η > = 0(2.24)

The ﬁve equations (2.1), (2.15), (2.17), (2.18), and (2.23) constitute the Kalmanﬁltering algorithm. The optimal weight K is called the Kalman gain matrix in Kalman13

ﬁltering. In optimal interpolation (OI), eqs.(2.1) and (2.15) are used with Pf constantin 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 deﬁned as
D

P =
'

δx(δx)>

En
1 X
δx(l) (δx(l) )>n − 1 l=1(2.25)(2.26)

where n and l denote the number of samples and a sample index respectively, and δx is avector of length N (the number of degrees of freedom of the model). An inﬁnite number ofsamples (n → ∞) 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 ẼU alsosatisﬁes eq.(2.27), where U is a unitary matrix, i.e., one that satisﬁes UU> = I. Eq.(2.27)
indicates that a given covariance can be perfectly represented by N samples
Ẽ = √hi1δx(1) · · · δx(N )N −1(2.28)

Usually N is O(107 ) in NWP models, which prohibits the explicit estimation of the covariance matrix. However, it is common for the covariance matrix P to be degenerate (seefor 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, andS is an N × m matrix whose columns are the corresponding eigenvectors. E is an N × m
matrix, which can be written as
E= √

h
i
1
δx(1) · · · δx(m)m−1(2.31)In this way, m samples express the covariance matrix P fairly well. As was the case forẼ, the choice of E is not unique because EU also satisﬁes eq.(2.30).
Using the square root of the covariance matrices, eq.(2.23) can be written as
Pfi= Efi (Efi )>= MPai−1 M>= 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 = M t .
Thus,
Efi = MEai−1(2.33)The TLM M gives the time evolution of the perturbation vector δx. By the deﬁnition of
the TLM, the original nonlinear model can approximate eq.(2.33), thus,MEai−1 =
'
'

i
h
1
a(1)
a(m)
Mδxi−1 · · · Mδxi−1
m−1
h
i
1
a(1)
a(m)
√
M (xai−1 + δxi−1 ) − M (xai−1 ) · · · M (xai−1 + δxi−1 ) − M (xai−1 )
m−1
i
h
1
a(m)
a(1)
√
(2.34)
M (xi−1 ) − x̄fi · · · M (xi−1 ) − x̄fi
m−1

√

15

where xa(l) and x̄f denote the lth ensemble member and the ensemble mean, respectively:
a(l)

xi

a(l)

= x̄ai + δxi(2.35)

x̄ai =

m
1 X
a(l)
x
m l=1 i(2.36)

x̄fi

m
1 X
a(l)
M (xi−1 )m l=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) canbe written
K = Pf H> (HPf H> + R)−1
= Ef (HEf )> [HEf (HEf )> + R]−1(2.38)

Thus, we can avoid storing the entire P matrix. We simply store Ef and HEf , 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

=
'
'

h
i
1
Hδxf (1) · · · Hδxf (m)
m−1
h
i
1
√
H(xf + δxf (1) ) − H(xf ) · · · H(xf + δxf (m) ) − H(xf )
m−1
h
i
1
√
H(xf (1) ) − H(x̄f ) · · · H(xf (m) ) − H(x̄f )(2.39)m−1√

Eq.(2.38) requires the inverse of a p × p matrix where p denotes the number ofobservations. If m < p, we can simplify computation of eq.(2.38) by expressing it as
K = Ef [I + (HEf )> R−1 HEf ]−1 (HEf )> R−1(2.40)

where the inverse of an m × m matrix is required. We assume here that diﬀerent observations are uncorrelated to each other, therefore the observational error covariance matrix
R is almost diagonal, and its inverse is trivial.
16

Finally, the Kalman ﬁltering 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 typesof ensemble update methods:
1. a perturbed observation (PO) method
2. a square root ﬁlter (SRF)
The PO method applies an independent data assimilation cycle to each ensemble member
yielding
Ea = (I − KH)Ef(2.41)

Pa = (I − KH)Pf (I − KH)>

(2.42)

This gives

where the analysis covariance is too small compared to eqs.(2.7) and (2.17). Thus, thisprocess 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 anadditional 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 = Ef T(2.43)

Then, eq.(2.17) yields
Pa = Ef TT> (Ef )>= (I − KH)Ef (Ef )>17(2.44)

By solving for T, we get an ensemble update formula. Since the choice of T is not unique(TU also satisﬁes eq.(2.44)), the ensemble update formula is not unique. Among theseveral 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 ﬁltering with SRF is realized by
the following procedures:
1. Ensemble forecast provides Ef and x̄f (eq.(2.34))
2. Generate Kalman gain matrix K using eq.(2.38) or (2.40)
3. Compute the analysis ensemble mean x̄a using eq.(2.1) with xa and xf replaced byx̄a and x̄f
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 ﬁltering (serial EnSRF)

The serial EnSRF (Whitaker and Hamill 2002) assumes an ensemble update of the form
Ea = (I − K̃H)Ef(2.45)In order for eq.(2.45) to give a solution of eq.(2.17),(I − K̃H)Pf (I − K̃H)> = (I − KH)Pf

18(2.46)

has to be satisﬁed (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 byhK̃ = Pf H> (HPf H> + R)−1/2

i> h(HPf H> + R)1/2 + R1/2

i−1

(2.47)(Andrews 1968), cf. eq.(10) in Whitaker and Hamill (2002). This formula allows forcorrelated observations. If observations are uncorrelated (R is diagonal) observations are
assimilated one at a time (serially), making the terms HPf H> and R scalar. In this case,eq.(2.47) can be simpliﬁed, and assuming K̃ = αK where α is a scalar value, one can get
the formula for α

Ã

α= 1+

s

R
f >
HP H + R!−1(2.48)

which was ﬁrst 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 required. In the present research, observations are assumed to be independent of eachother, which makes only the computation of eq.(2.48) necessary. Thus, K and H are vectors 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 matrix P, where m is O(102 ). This approximation is not usually justiﬁed
because the P reproduced from the m members of ensemble contains large sampling errors,
resulting in artiﬁcially large covariances between distant points. Thus, it is necessary to
localize the inﬂuence of observations. Intuitively, covariances between distant points canbe assumed to be zero, since it is reasonable to assume there is no correlation between, forexample, 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 thecenter point. In this way, the input function f is localized around the center point. The
factor S(r) is given by a ﬁfth-order piecewise rational function:
1 − 14 r5 + 12 r4 + 58 r3 − 53 r2






S(r) =(r ≤ 1)

1 5
1 4
5 3
5 2
2 −1

12 r − 2 r + 8 r + 3 r − 5r + 4 − 3 r






 0

(1 < r ≤ 2)

(2.49)(2 < r)cf. eq.(14) in Hamill et al. (2001). Here, r is deﬁned as
d
r=p
10/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 distantpoints. In serial EnSRF, the Kalman gain K, now an N dimensional vector, is localizedaround 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 arecaused by sampling errors. This method of ”gain localization” is equivalent to covariancelocalization (cf. eq.(5) of Houtekamer and Mitchell 2001) as follows³´−1

K = (S ◦ Pf )H> H(S ◦ Pf )H> + R³≈ S ◦ Pf H> (HPf H> + R)−1

(2.51)´(2.52)

Here, we used the fact that HPf H is just a scalar in the serial treatment of observations.
If the observation is at a grid point, eq.(2.52) becomes exact.
20

Gaussian−like weighting function (sigma=10)
1
Gaussian
Gaussian−like weighting

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
−50

−40

−30

−20

−10

0

10

20

30

40

50

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 ﬁltering (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 distantpoints. 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 = GP̂f G>(2.53)

deﬁnes a variable transformation G. P̂f is a diagonal matrix with just r components
diﬀerent from 0. Ignoring degenerate eigenvectors, we can deﬁne G and P̂f 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 isgiven byx = Gx̂(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 byx̂a = x̂f + P̂a Ĥ> R−1 d
P̂a =

h

(P̂f )−1 + Ĥ> R−1 Ĥh= P̂f I + Ĥ> R−1 ĤP̂f(2.57)
i−1(2.58)

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 byÊa = Êf Y(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
P̂a = Êf YY> (Êf )>(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 byh

Y = I + (Êf )> (P̂f )−1 (P̂a − P̂f )(P̂f )−1 Êf

i1/2(2.62)

cf. eq.(41) in Ott et al. (2004). This choice of Y yields the minimum diﬀerence between
the forecast ensemble perturbations and the analysis ensemble perturbations. This Y isadvantageous to the smoothness and physical balance of the analysis ﬁelds.

2.2.6 Observational error covariance localization
In the LEKF formulation described in the previous section, LEKF localizes around ananalyzed point with Heaviside-like step weighting whereas serial EnSRF used Gaussianlike weighting function. Hunt (pers. comm.) suggested a Gaussian-type error covariance23

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 tocarry out in the LEKF formulation) by weighting observational error according to the
distance from the center of the local patch. Intuitively,³´−1

K = (ρ ◦ Pf )H> H(ρ ◦ Pf )H> + R³´−1' Pf H> HPf H> + (ρ−1 ◦ R)

(2.63)(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 componentsof 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 weapply 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, respectively. We call this localization ”observational error covariance localization” originallysuggested by Hunt (pers. comm.).

24

2.2.7

Covariance inﬂation

Kalman ﬁltering 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 ﬁltering, so thatobservations are given enough weight. If not, the process ampliﬁes feedback and eventually
results in ﬁlter divergence.
A simple way to inﬂate covariance is to multiply the covariance matrix Pf or Pa bya factor slightly larger than 1
Pf ← Pf × (1 + δ)

(2.66)

where δ is a small positive number (δ << 1). Using the inﬂated covariance Pf , observations have larger eﬀects and we prevent ﬁlter divergence. This operation is called”multiplicative covariance inﬂation”.
Ott et al. (2004) proposed a more sophisticated way to inﬂate covariance in thecontext of LEKF, which is called ”enhanced variance inﬂation”. The diagonal componentsof the analysis error covariance matrix in the ”hat space” P̂a are inﬂated as follows:P̂a ← P̂a +δ × trace(P̂a )I
r(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¶µ1
P robf (x) ∝ exp − (x − xf )> B−1 (x − xf )2(2.68)In the same way, the PDF with respect to the observation is given byµ¶1
P robo (x) ∝ exp − (Hx − yo )> R−1 (Hx − yo )2(2.69)The joint probability is written as
P robb&o (x) = P robb (x) · P robo (x)
µ

¶

1
1
∝ exp − (x − xb )> B−1 (x − xb ) − (Hx − yo )> R−1 (Hx − yo )22(2.70)

Since data assimilation is a maximum likelihood estimate, the analysis state x is obtainedby maximizing eq.(2.70). This is equivalent to minimizing the cost function

J(x) =

1
1(x − xb )> B−1 (x − xb ) + (Hx − yo )> R−1 (Hx − yo )22

(2.71)

At a minimizer x∗ , the gradient of the cost function is equal to 0 , that is,∇J(x∗ ) = 0(2.72)

3DVAR algorithms look for the solution that satisﬁes 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 theforecast and observations.
In the incremental form, the cost function (eq.(2.71)) can be rewritten as
1
1
J(xb + δx) ≈ δx> B−1 δx + (Hδx − d)> R−1 (Hδx − d)22(2.73)

where we deﬁned 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 thegenerally 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∇J(δx) = B−1 δx + H> R−1 (Hδx − d) = 0(2.76)

Solving eq.(2.76), we getδx = (B−1 + H> R−1 H)−1 H> R−1 d(2.77)= BH> (R + HBH> )−1 d(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 withlinearized H. Since 3DVAR uses the time-independent background error covariance B obtained 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 deﬁned asdxi= xi−1 (xi+1 − xi−2 ) − xi + F
dt(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 ﬁrst 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 diﬀerencing 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 thesame problem with diﬀerent 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 sameform, and the optimal weight (eq.(2.15)) is computed explicitly. To obtain the backgrounderror covariance B, we apply the NMC method (Parrish and Derber 1992), accumulating28

diﬀerences between 24-hour forecasts and 18-hour forecasts valid at the same time as
samples of background errors.
In this simple system, the ﬁve equations for Kalman ﬁltering can be solved explicitly. Thus, the ”full Kalman ﬁltering” (full-KF) without low-rank assumption has been
implemented.
As for the serial EnSRF and LEKF, model-independent core modules have beendeveloped. 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 incore modules is set to 40. The core subroutines are called in the loop of observations
because observations are assimilated serially. After assimilating each observation, thebackground ensemble is updated serially. After assimilating all observations, we retain theﬁnal 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

l

C

l2

Figure 2.2: Deﬁnition of the local regions around a center point (C) for LEKF in theLorenz-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 ﬁrst 60 days in the analysis
experiments are considered as the spin-up period. To obtain the 1-year (365-day) meanof 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 deﬁnes 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 ﬁltering methods show the same performancein 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 parametersettings 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 ﬁxed at 10.

twice as large (RMSE=0.40).
When serial EnSRF uses 40 ensemble members, the ﬁlter is stable even without
the Schur product to localize the Kalman gain. When ensemble size is reduced, theSchur 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 isshown graphically in Fig.2.3, which is similar to Fig.3(b) of Whitaker and Hamill (2002).
The localization length scale σ is diﬀerent from the covariance ﬁlter length scale L in
Whitaker and Hamill (2002). Whitaker and Hamill (2002) deﬁned the length scale L so
the covariance becomes 0 when d > L where d denotes distance. Our deﬁnition is given
by eqs.(2.49) and (2.50), where the distance d that yields zero covariance is given by
r

d>2×

10
×σ
3

(2.80)

Thus, the relationship between the length scale L in Whitaker and Hamill (2002) and thelength scale σ used in this experiments is given asr

L = 2×' 3.65σ31

10×σ
3(2.81)

Lorenz96 Analysis RMSE NOBS=40 Serial EnKF (NBV=10)
0.1

5
0.2
27
0.

0.325

0.25

0.07

0.

0.2

5

22

5

0.225

0.05

0.27

5

0.06

0.3

covariance inflation factor

0.08

5

0.3

0.09

25

2
0.

0.04

0.21

0.02

FILTER DIVERGENCE

0.21

25

0.2

0.3

0.275

0.25

0.03

0.01

1

2

3

4

5
6
Schur product length scale

7

8

9

10

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 lengthscale σ and the covariance inﬂation paramter δ, respectively. The minimum error 0.20 isobserved when σ = 5 and δ = 0.02. ”FILTER DIVERGENCE” denotes the region with
RMSE of more than 1.0.

In considering the diﬀerent measure of the length scale, Fig.2.3 shows almost identical
performance as Fig.3(b) in Whitaker and Hamill (2002), which supports the correctnessof 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 ﬁgure as Fig.2.3 in the case of LEKF. LEKF has three parameters: 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Lorenz96 Analysis RMSE NOBS=40 LEKF (NBV=10)

0.3

0.
27
5

0.1

0.09

0.25

0.25

0.07

25

0.06

5

0.27
5

0.2

0.05

25

0.2

0.

covariance inflation factor

0.08

25

0.04

0.275

0.03

2
0.

0.01

0.21

0.25

0.02

1

5
22

1
0.2

0.
2

3

4

5
6
Size of local patch (l)FILTER DIVERGENCE
7

8

9

10

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 ﬁxed to 1.
The minimum error 0.20 is observed when l = 5 and δ = 0.01.

has been observed with respect to l2 , l2 is ﬁxed to 1. Note that the direct comparison
between serial EnSRF and LEKF may be misinterpreted since the way of localization is
diﬀerent between serial EnSRF and LEKF, so that the metric in the horizontal axis may
be diﬀerent. 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 ﬁlter (i.e., the farther
away from the region of ﬁlter divergence).
• The larger the covariance inﬂation factor, the more stable the ﬁlter.

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 ﬁltering.
• 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 ﬁlter 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 ﬁxed 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 backgrounderror 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 islarger than the observational error, whereas those of Kalman ﬁlters are much smaller. The
diﬀerence between 3DVAR and Kalman ﬁlters becomes larger here, thus, Kalman ﬁltersshow a larger advantage with fewer observations. All the Kalman ﬁltering methods show
the same performances (RMSE=0.33) with their optimal settings.
Fig.2.5 shows the same ﬁgure as Fig.2.3 but in the case of 20 observations. Comparedto the case of 40 observations, the ﬁlter becomes considerably more unstable and the region”FILTER DIVERGENCE” becomes larger, although the basic characteristic is similar.
34

Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=10)

0.4

5

45

42

0.

0.

0.4

75

0.1

0.09

0.5

0.3

75

0.3

75

0.08

4

0.07

0.05

0.425

0.475
0.45

0.35

0.5

0.06

5

37

0.

0.375

0.04

0.35

0.4

0.03

FILTER DIVERGENCE

0.45
0.475

0.425

covariance inflation factor

0.

0.02

0.5
0.01

1

2

3

4

5
6
Schur product length scale

7

8

9

10

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 ﬁgure as Fig.2.5 but ensemble size is increased to 20. Whenensemble size is increased, the ﬁlter becomes more stable. In addition, the error becomesslightly 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 thesame ﬁgure as Fig.2.5 but with 8 ensemble members. Fig.2.7 shows a similar structure as
Fig.2.5 but with larger ”ﬁlter divergence” area and with larger errors almost everywhere.
Smaller ensemble size requires smaller covariance length scales. The minimum error 0.34 isobserved when σ = 5 and δ = 0.04, which is slightly worse than in the case of 10 ensemble
members.

35

Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=20)
0.
5

0.1
45
0.

0.09

0.425

25

0.475

0.4

0.4

5
0.37

0.4

0.07

0.05

75

0.35

0.35

5

0.32
0.3

0.45

0.375

0.03

0.01

1

25

FILTER DIVERGENCE
0.

32

5

0.4
5
0.42

0.02

0.3

5

0.04

0.3
0.4

0.425

0.06

0.45

0.5

0.375

0.475

covariance inflation factor

0.08

2

3

4

5
6
Schur product length scale

7

8

9

10

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

Lorenz96 Analysis RMSE NOBS=20 Serial EnKF (NBV=8)
0.1

0.09

5
0.4

0.5
0.475

4

0.

5
42
0.

0.4

5

0.37

5
0.37

0.3

5

0.3

5

0.05

0.45

0.06

0.4

25

0.4

0.07

0.475

covariance inflation factor

0.08

0.04

FILTER DIVERGENCE

0.4

0.01

0.45
0.475

0.02

0.42
5

0.03

1

2

3

4

5
6
Schur product length scale

7

8

9

10

Figure 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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 ﬁxed as l = 6. The other two parameters, l2 and δ, are changed.

Fig.2.8 shows the same ﬁgure in the case of LEKF. Here, l2 is ﬁxed to 0. It seems
LEKF has a smaller ”ﬁlter 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 ﬁxed 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 theruns may not be suﬃcient to get reliable RMSE values, although the same computations
with diﬀerent random numbers for observational errors performed several times yielded
the same overall results. Table 2.3 suggests that the LEKF errors are rather insensitive tol2 , 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 ﬁlter becomes divergent (Fig.2.8).

38

Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8)
0.1
0.5

5

5

0.4
75

0.45

0

0.5
2

0.09

0.42
5
.42

4

0.4

5

45

7
0.3

0.

0.05

0.04

FILTER DIVERGENCE
0.42
5

0.03

5 0.5
0.01

1

0.45

7
0.4

0.02

0.

0.42
5

0.5
0.4

75

0.06

0.425

0.07

0.4

covariance inflation factor

0.

4

0.08

2

3

4

5
6
Size of local patch (l)

7

8

9

10

Figure 2.8: The same as Fig.2.4 but LEKF with 8 ensemble members and enhanced
variance inﬂation. l2 is ﬁxed to 0. The minimum error 0.35 is observed when l = 5 andδ = 0.03.

39

Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8)

0.

42

0.5

0.47

5

5

0.425

0.08

0.45

0.07

0.05

0.42
5

0.4
75

0.
5

0.06

0.425

covariance inflation factor

5

0.55

0.09

0.4

0.475
45
0.

0.525

0.1

0.4

0.04
FILTER DIVERGENCE
0.03

0.4
0.

0.02

0.01

1

5

0.425
5

0.47
5

2

3

4

5
6
Size of local patch (l)

7

8

9

10

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 ﬁgure as Fig.2.8 but without enhanced variance inﬂation,
replaced by the same multiplicative covariance inﬂation used for the serial EnSRF. Enhanced variance inﬂation contributes to stabilize the ﬁlter, the minimum error becomes
smaller with the smaller inﬂation factor. In addition, with enhanced variance inﬂation, theerrors are smaller than with multiplicative inﬂation almost everywhere. Thus, enhanced
variance inﬂation contributes positively to the ﬁlter performance.
We also applied observational error covariance localization within LEKF. Fig.2.10
shows the same ﬁgure as Fig.2.9 but with observational error covariance localization.
Here, we ﬁxed local patch parameters to l = 15 and l2 = 0, and varied localization scale
σ of the observational error covariance localization. With observational error covariance

40

0.5
0.475
0.4
5

0.375

0.4

0.09

0 .4
25

Lorenz96 Analysis RMSE NOBS=20 LEKF (NBV=8) w/LOCALIZATION
0.4

0.1

37
0.

0.3

75

0.07

0.

35

5
0.3

0.4

0.05

0.5
0.475
0.45

0.06

0.425

covariance inflation factor

5

0.08

0.04

0.375

FILTER DIVERGENCE

0.03
0.425
5
0.4

0.02

0.3

5

4

0.

0.01

1

2

3

4

5
6
Schur product length scale

7

8

9

10

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 ﬁxed tol = 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 ofenhanced 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 thepractical 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:
DE

dd> = (1 + δ)HPf H> + R(2.82)

a similar equation is found in eq.(3) of Houtekamer et al. (2005). Here, d denotes theobservational increment deﬁned 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 ﬁltering in a model output statistics (MOS) (see for example,
appendix C of Kalnay 2003), we can apply Kalman ﬁltering to make an online estimationof δ. 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:δbσo + δoσbσb + σo
!
Ã
σb
σb
=
1− b
σ + σoδa =(2.84)

σa(2.85)

This analyzed δ a can be used for the variance inﬂation parameter at the current step. δ aand σ 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 both42

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 isapplied. l2 is ﬁxed to l2 = 0. In the cases with observational error covariance localization,l is ﬁxed 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 observedwith 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 covarianceinﬂation adaptively: In a global model the use of the same inﬂation for both tropical andmidlatitude 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 isshown in Table 2.5. The computation is a whole 425-day forecast-analysis cycle, containing
1700 data assimilations. The only diﬀerence between 3DVAR and full-KF is covariance43

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 isslower 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 independent, 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 local ensemble transform KF (LETKF, Hunt 2005) avoids the eigenvalue decomposition andis 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 diﬀerent 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 havebeen able to reproduce the results of Whitaker and Hamill (2002) and of Ott et al(2004).
• All the Kalman ﬁltering methods (full-KF, serial EnSRF, and LEKF) outperform
44

3DVAR, and the advantage becomes larger with fewer observations.
• All the Kalman ﬁltering methods (full-KF, serial EnSRF, and LEKF) have comparable performance with their optimal settings.
• Observational error covariance localization contributes to better LEKF performancewith little additional computational cost.
• Enhanced variance inﬂation in LEKF contributes to better ﬁlter 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 thanserial EnSRF with similar parameters.
• Smaller localization length scales and larger covariance inﬂation factors stabilize
EnKF, although both give suboptimal ﬁltering.
• 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 constantinﬂ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 fullKF) 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. What is the relative advantages and disadvantages between serial EnSRF and LEKF?
3. How sensitive is the data assimilation to experimental settings?
4. What are the characteristics of the analysis and forecast error ﬁelds?
A goal of this chapter is to investigate and compare the performance of the threemethods. In addition, we perform sensitivity experiments with diﬀerent 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 ﬁlter performance. We also see the characteristics of the46

analysis and forecast error ﬁelds, suggesting how EnKF works.

3.1.2

The SPEEDY model

The SPEEDY model (Molteni 2003) is a recently developed atmospheric general circulation model (AGCM) with a spectral primitive-equation dynamic core and a set of simpliﬁed physical parameterization schemes (SPEEDY stands for Simpliﬁed Parameterizations, primitivE-Equation DYnamics). The goal of this model is to achieve computationaleﬃciency while maintaining characteristics similar to the state-of-the-art AGCMs withcomplex 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 magnitudeless 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 atmosphericcirculation 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 usedin more complex GCMs, such as convection (a simpliﬁed mass-ﬂux scheme), large-scalecondensation, clouds, short-wave radiation (two spectral bands), long wave radiation (fourspectral bands), surface ﬂuxes of momentum and energy (bulk aerodynamic formula),
and vertical diﬀusion. Details of the simpliﬁed physical parameterization schemes of theSPEEDY model can be found in Molteni (2003), especially in its Appendix which isavailable 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 speciﬁed by monthly means, and bare-surface albedo
and fraction of land-surface covered by vegetation, which are speciﬁed by annual-meanﬁelds. 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 ), speciﬁc 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 insigma levels, but outputs are in both sigma levels and pressure levels. The speciﬁc heightsof 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 experiments, 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 foreach 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 asmodel 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 impossible to implement the variational formulation explicitly because of the large number ofdegrees of freedom of the system. Usually, NWP models have prognostic variables withdimensions 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 ﬁfth-generation Pennsylvania State University/National Center for
Atmospheric Research Mesoscale Model (MM5). We follow a similar approach as Barker
et al.
We deﬁne a variable transformation U as follows:δx = Uδv(3.1)

Here, δx is deﬁned by eq.(2.75), δv is a vector with the same dimension as δx. Substituting
eq.(3.1) into eq.(2.73), we getJ(δv) =

1 > > −11
δv U B Uδv + (HUδv − d)> R−1 (HUδv − d)
2
2

(3.2)

In order for the covariance matrix to be the identity,
B = UU>(3.3)

should be satisﬁed. Then, eq.(3.2) can be rewritten as

J(δv) =11 >δv δv + (HUδv − d)> R−1 (HUδv − d)
2
2

(3.4)

Diﬀerencing eq.(3.4) with respect to δv, we get the gradient of the cost function∇J(δ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)) andits gradient (eq.(3.5)), a quasi-Newton minimizer ﬁnds 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 algorithm. 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 correlationand the inter-variable error correlation, respectively. Moreover, it is assumed C can beseparated into horizontal and vertical correlations. Usually, the horizontal correlation isassumed 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 samevertical 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 independent of each other. Some prognostic variables are strongly dependent on each othermostly 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 thesystem.

3.2.2

Implementation on the SPEEDY model

In the previous section, we described a general idea to implement 3DVAR on a practicalatmospheric model following the approach by Barker et al. (2004). Here, we describe thedesign of the variable transformation U in the present implementation.
The ﬁrst part is the error standard deviation A in 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 requiredby considering the full spatial dependence.
For the spatial error correlation C in eq.(3.6), the spatial correlation is separated into
vertical and horizontal correlations. The vertical resolution of the SPEEDY model is socoarse 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 ﬁlter (RF) technique (Purser et al. 2003) has
been implemented to estimate the Gaussian correlation shape. The 4th order RF is aquasi-Gaussian ﬁlter composed of a pair of forwarding and backwarding processes
Bi = βAi + α1 Bi−1 + α2 Bi−2 + α3 Bi−3 + α4 Bi−4

(3.7)

Ci = βBi + α1 Ci+1 + α2 Ci+2 + α3 Ci+3 + α4 Ci+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 ﬁltered value
B is obtained. In eq.(3.8), that is the backwarding process in the descending order of i, B
is input and the ﬁltered 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 β andα can be obtained by solving a 4th degree equation according to the correlation lengthscale. The details of the recursive ﬁlter 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 − r1 ug (ps , T )(3.9)

vu = v − r2 vg (ps , T )(3.10)

where r1 and r2 are regression coeﬃcients determined from statistics, and ug and vg are52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)):
f k × vg = −RT ∇ ln ps − ∇φ

(3.11)

where f , k, R, and φ denote Coriolis parameter, vertical unit vector, gas constant, andgeopotential height, respectively. Other variables (T , q, and ps ) are assumed to be independent of the unbalanced wind (uu and vu ). Thus, only the geostrophic balance isconsidered for inter-variable correlation.

3.2.3

Background error statistics - NMC method

The design of variable transformation U described in the previous section requires background error statistics. We need to specify the error standard deviation, the length scaleof horizontal error correlations, and the regression coeﬃcients from the statistics. Background error statistics were computed using the NMC method (Parrish and Derber 1992),
accumulating the diﬀerence between 24-hour forecast and 18-hour forecast as a sample ofthe background error. In the present experiment, 87 samples from January 10 to January
31 were accumulated for the statistics. The analysis ﬁelds for the ﬁrst time are obtainedusing 3DVAR with a reasonable guess of background error statistics, such as 500km for thehorizontal 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 backgrounderror 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 uwind (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, areobserved in mid-latitudes.

geostrophy is observed in the mid-latitudes especially at the jet levels. The regression coefﬁcients 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 spatial dependence, especially in zonal structures in error standard deviation. Fig.3.3 shows
the vertical structure and latitudinal dependence of horizontal error correlations. Thereis 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 physicallymore dense in higher latitudes.
The ﬁgures shown as background error statistics so far are only small parts of the54

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 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 showcorrelation.

56

background error information that are considered in the present 3DVAR. The present
3DVAR considers all components of standard deviation, the standard deviation has thesame 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 latitudeand 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 tovary 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 verticalcorrelations averaged horizontally. Fig.3.4 shows vertical correlations using the NMC
method. Some vertical correlations are seen at lower levels, but almost no correlation isobserved at upper levels. Except for moisture, only the ﬁrst and second levels have large
correlations (about 0.4), other levels show less than 0.2 correlations. Thus, the verticalindependence is not a bad assumption in 3DVAR. Moisture, shown by long-short dashedlines with black square, has largest vertical correlation at levels 2, 3, and 4, which may berelated 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 incrementof 1.0m/s. Fig.3.5 shows responses when the observational station is located at midlatitude (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) when57

Figure 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

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 istested using the dense observational network (Fig.3.9). Fig.3.7 shows analysis increment(top panel) and the diﬀerence between truth and ﬁrst guess (bottom panel). Analysis
increment has a very similar structure as the diﬀerence between truth and ﬁrst 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 considered. The boundary is cyclic in longitudinal direction, whereas in latitudinal direction, theboundary 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 diﬀerent 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 theobservational 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 withthe observational increment of 1.0m/s at mid-latitude at the 4th level. Analysis incrementof 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 diﬀerence between truth
and ﬁrst guess (bottom panel). The ﬁelds are u-wind ﬁelds (in m/s) at the 4th level.

62

l

l2

Figure 3.8: Deﬁnition 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 fromthe 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 inthe case of the Lorenz-96 model. The region is a square in the grid space, the physicalshape is not a square because of the skewness of the map in higher latitudes. Since thelocal patches are overlapped, only the values inside the smaller region denoted by l2 areaveraged to obtain the global analysis.
63

Cyclic boundary conditions are applied in the longitudinal direction. In the latitudinal 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 diﬀerent 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 aﬀect other variables
at all levels. However, only the variables at the bottom level aﬀect ps . In this way, thecombination 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 developed so that it can consider vertical correlations. For comparison with 3DVAR andserial EnSRF, vertical correlation is not considered, but the impact of vertical correlationis shown later.

3.4 Experimental setup
So far, we described the implementation of three methods (3DVAR, serial EnSRF, andLEKF). 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]

Speciﬁc 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 thesparse 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 inWestern 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 nointerpolation 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 perfectly as shown later. Thus, in the dense network, there would be no signiﬁcant diﬀerence
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 percentsof 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 percents 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 expected 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 theinitial date for data assimilation experiments after 1-year integration. We use the same67

nature run for all data assimilation experiments in this chapter.
We generate another independent nature run starting from the atmosphere at rest
at a diﬀerent time. The diﬀerence between the two nature runs temporally averaged in
February, 1982 is 106.37 meters measured by the RMS diﬀerence of 500hPa height ﬁeld,
the error magnitude without data assimilation. We choose ﬁelds at arbitrary times fromthe second nature run to obtain the initial conditions of data assimilation. For 3DVAR,
a ﬁeld is chosen. For EnKF, we choose ﬁelds at arbitrary times for the initial ensemble.
The same initial ensemble is used for serial EnSRF and LEKF. If the ensemble size is
inﬁnitely large, the ensemble mean becomes the climatology. Thus, this is equivalent tobeginning the data assimilation cycle with climatology in 3DVAR. This choice is practicalin 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 ﬁeld (”FORECAST”) by the SPEEDY model gives the ﬁrst guess ﬁeld for thedata assimilation system (”DA”), and the true ﬁeld (”TRUTH”) with random noise provides observations. ”DA” may be 3DVAR, serial EnSRF, or LEKF. The output of ”DA”
is the analysis ﬁeld (”ANALYSIS”) that provides the initial condition for the next model
integration. ”ANALYSIS” contains high frequency components, whereas ”INITIAL”, theoutput by the SPEEDY model at forecast time 0, does not because the spectral transformation of the model smoothes out the high frequency components. Since the SPEEDY
model does not see the high frequency components, the smoothed ”INITIAL” ﬁeld is considered as the analysis ﬁeld for veriﬁcation. Thus, the analysis error is deﬁned as thediﬀerence between ”INITIAL” and ”TRUTH”.
We use the multiplicative covariance inﬂation (eq.(2.66)) both for serial EnSRF and

68

FORECAST

SPEEDYINITIAL

Verification

ANALYSIS

DA

TRUTH

Add random noise as observational errors

Figure 3.12: Schematic of data assimilation experiments on the SPEEDY model. ”DA”
denotes the data assimilation system. The ”initial” ﬁeld is the output by the SPEEDY
model which is diﬀerent from the ”analysis” ﬁeld because the high frequency components
are smoothed out by the spectral transformation of the model. The diﬀerence between”initial” and ”truth” is deﬁned as the analysis error.

69

LEKF. The enhanced variance inﬂation (eq.(2.67)) has been tested, but the ﬁlter quickly
diverges in the SPEEDY model. The inﬂation factor is ﬁxed 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 andensemble 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 ﬁeld. The solid line shows the resultin the case of a dense observational network. The dense observational network is denseenough for 3DVAR to perform almost perfectly; the analysis RMSE of 500hPa height isless than 5 meters, much smaller than the operational NWP system for the real atmosphere. 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, theperformance 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 regeneratedusing the sparse network. The new background error information for the sparse network

70

Figure 3.13: Analysis RMSE of 500hPa height ﬁeld in the case of 3DVAR with the dense
observational network (solid line), sparse observational network using the same backgrounderror 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 (dotdashed line).

gives worse 3DVAR performance than the original background error information basedon the dense network (dot-dashed line). The three lines of 3DVAR with the backgroundstatistics 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
deviation in space. Fig.3.14 shows analysis RMSE of 500hPa height in the case of the sparse
observational network with diﬀerent 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 dashedline shows the results in the case of horizontal mean of standard deviation, in which casethe standard deviation is a function of only the vertical levels and variables. The dottedline shows the results in the case of zonal mean of standard deviation, in which case, thestandard 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 thatthe full spatial structure of the error covariance obtained with 87 samples is dominatedby sampling errors.
Since we observed the best performance (31.14 meters) when we multiplied thebackground error standard deviation by 2.0 and took the horizontal mean of the originalerror variance obtained from NMC-method statistics, hereafter, we use this conﬁguration
as an optimal setting of 3DVAR.
In order to see the error distributions in space, analysis RMSE of 500hPa height ﬁeld
is plotted in space. Analysis errors are averaged only in time to obtain the spatial distribution. Fig.3.15 shows the horizontal structure of analysis RMSE of 500hPa height ﬁeld.
Shades and contours show analysis RMSE and mean analysis ﬁelds, 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 ﬁnd the lattice-like pattern. We

72

Figure 3.14: Analysis RMSE of 500hPa height ﬁeld 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 withzonally uniform background error standard deviation (dotted line).

RMSE[m]3DVAR(full)

3DVAR(zonal mean)

3DVAR(hmean)

NO ASSIM

36.65

34.62

31.14

106.37

Table 3.3: Analysis RMSE of 500hPa height ﬁeld 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 initial one-month spin-up period. ”3DVAR(full)”denotes the case that full spatial dependence of background error standard deviation isconsidered. ”3DVAR(zonal mean)” and ”3DVAR(hmean)” denote the cases that zonallyand 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 ﬁeld in the case of 3DVARwith horizontally uniform background error standard deviation and the sparse observational network. Shades and contours show analysis RMSE (in meters) and mean analysis
ﬁelds (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 ﬁeld 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 regularbias structure. If we subtract error bias from the analysis error ﬁeld and recalculate theanalysis 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 ﬁeld. Although
this ﬁeld 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 ﬁeld (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 ﬁelds. 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 ﬁelds, respectively.
If we ignore the large errors in the polar regions, we see large analysis error in the tropicalregions around the 200hPa level in both ﬁgures. In addition, both ﬁgures show striped
pattern, where the minima correpond to observational locations. For the temperatureﬁeld (bottom panel), we see the tendency for larger errors to appear in upper levels in thetropics. 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 networkshows almost perfect performance: 3.58m RMSE of 500hPa height ﬁeld.
• The 3DVAR data assimilation cycle experiments using the sparse observational network show much worse performance: 31.14m RMSE of 500hPa height ﬁeld.
• 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 agood 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 tothe observational locations.
76

Figure 3.17: Zonal structure of analysis RMSE of u-wind ﬁeld (top panel, in m/s) andtemperature ﬁeld (bottom panel, in K) in the case of 3DVAR with horizontally uniformbackground error standard deviation and the sparse observational network. Shades andcontours show analysis RMSE and mean analysis ﬁelds, 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 ﬁeld in meters using serial EnSRF on theSPEEDY 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 ischosen as integers from 1 to 6. ”D” denotes ﬁlter divergence. 4% multiplicative covarianceinﬂation is applied.

3.5.2

Serial EnSRF

Table 3.4 shows the analysis RMSE of 500hPa height ﬁeld. The ensemble size is chosen
to be 10, 20, and 30, and the horizontal correlation length scale for the Schur product ischosen 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 thehorizontal scale of 3.0, which shows almost perfect performance (as good as 3DVAR withthe dense observational network).
A smaller ensemble size requires a smaller horizontal length scale of the Schur product. 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 ﬁeld in78

Figure 3.18: Analysis RMSE of 500hPa height ﬁeld 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.

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, thatgives 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ﬁeld 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 ﬁeld, 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 ﬁeld, similar to Fig.3.15, inthe case of serial EnSRF with 20 ensemble members and localization scale of 2.0. Shades
show analysis RMSE (in meters), contours show mean analysis ﬁelds (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 withlow wavenumbers exist in polar regions. We ﬁnd 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 ﬁeld 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 ﬁeld (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 ﬁelds.
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 ﬁelds. 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 ﬁelds, 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 ﬁelds. 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 solidline) at all levels. Fig.3.21(d) shows the analysis RMSE of speciﬁc humidity ﬁelds. The
order of the performances is basically the same as the previous ﬁgures (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 diﬀerence among diﬀerent data assimilation
schemes and settings in higher altitudes. Lower altitudes contain more humidity, andmostly the RMSE exceeds the observational error standard deviation (10−4 kg/kg, shownby the thin solid line).
Fig.3.22 shows the zonal structures of analysis RMSE of the u-wind (top panel) andtemperature (bottom panel), corresponding to Fig.3.17 of 3DVAR. Shades and contours
show analysis RMSE and mean analysis ﬁelds, 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 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 speciﬁc humidity (q) ﬁelds, 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 ﬁeld
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 theLorenz-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 latticelike pattern to the analysis error bias.
• Not only 500hPa height ﬁelds, but also other ﬁelds at other levels show much better
performance than 3DVAR.

3.5.3

LEKF

Table 3.5 shows the analysis RMSE of 500hPa height ﬁeld, similar to Table 3.4 in the caseof 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 ﬁeld (top panel, in m/s) andtemperature ﬁeld (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 andmean analysis ﬁelds, 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 ﬁeld in meters using LEKF 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 local patch parameter l is chosen as integers from 1 to 6.
”D” denotes ﬁlter 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 theanalysis RMSE of 500hPa height ﬁeld in the cases of LEKF with 10 (solid line), 20 (dashedline), 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ﬁeld is plotted in space, similar to Figs.3.15, 3.16, 3.19 and 3.20. Fig.3.24 shows thehorizontal structure of analysis RMSE of 500hPa height ﬁeld, similar to Figs.3.15 and86

Figure 3.23: Analysis RMSE of 500hPa height ﬁeld 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, theRMSE of 3DVAR is shown by the short dashed line.

873.19. Large errors exist in polar regions where no observations are available, but theerrors are slightly larger than those in serial EnSRF. We can see a lattice-like patterncorresponding 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 thesame 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 ﬁeld shows almost no latticelike 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 ﬁgures. Thus,
3DVAR, LEKF with 20 and 30 ensemble members are shown by short dashed lines, solidlines, 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 ﬁelds. Smaller RMSE is observed at thelevels 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, shownby 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 ﬁelds,
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 ﬁeld, 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 ﬁelds (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 ﬁeld (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
ﬁelds. Dots indicate observational locations.

90

there is almost no diﬀerence 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 ﬁelds. Generally, higher altitude gives smaller
RMSE, that is clear in 3DVAR. Almost all cases of LEKF provide smaller RMSE thanobservational 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 speciﬁc humidity ﬁelds,
where we can see that higher altitude is drier, so the RMSE is smaller in higher altitudes.
All show similar performances, and no clear diﬀerence can be seen, whereas there was clear
diﬀerence in the case of serial EnSRF (cf. Fig.3.21(d)). Lower altitudes contains morehumidity, and the RMSEs exceed the observational error standard deviation (10−4 kg/kg,
shown in the thin solid line).
Fig.3.27 shows the zonal structures of analysis RMSE of the u-wind (top panel) andtemperature (bottom panel), similar to Figs.3.17 and 3.22 of 3DVAR and serial EnSRF.
Shades and contours show analysis RMSE and mean analysis ﬁelds, respectively. We cansee the same basic features seen in 3DVAR and serial EnSRF: large errors in u-wind ﬁeld
exist in upper tropics, large errors in the temperature ﬁeld 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 createdin the same way as Table 3.5, that is, the temporal mean of the analysis RMSE of 500hPa
height ﬁeld is computed using various parameter settings. Here, the local patch parameter
is ﬁxed to l = 2, and the ensemble size is ﬁxed 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 speciﬁc
humidity (q) ﬁelds, respectively. The observational error standard deviations are shownas thin solid lines wherever applicable. The localization length scale of the Schur productis 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 ﬁeld (top panel, in m/s) andtemperature ﬁeld (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 ﬁelds, 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 ﬁeld 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 ﬁxed
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 ﬁeld 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. However, 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 ﬁelds, but also most other ﬁelds 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 isabout 35, which requires about 5 seconds of computational time.
Table 3.7 shows timing of serial EnSRF for each analysis. Here, the sparse observational 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 observations 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 thenumber of observations. The system has been accelerated by optimizing the code as much

95

NBV

10

20

Celeron 2.7GHz (Northwood)

2:20

4:30

Celeron D 2.53GHz (Prescott)

1:15

2:20

30

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)

Celeron 2.7GHz (Northwood)

2:50

8:30

Celeron D 2.53GHz (Prescott)

1:35

6:00

30(4)

16:30

Table 3.8: Timing (min:sec) of LEKF with 10, 20, and 30 ensemble members on theSPEEDY model when the sparse observational network is applied. The local patch parameter l is chosen as an optimal setting shown in the parenthesis, and l2 is ﬁxed to0.

as possible. However, we do not apply the new technique to avoid eigenvalue decomposition(local ensemble transform Kalman ﬁlter, LETKF) proposed by Hunt (2005). Note that
LEKF can be parallelized very eﬃciently, thus, LEKF can be accelerated almost perfectlyon parallel computers.
Table 3.8 shows timing of LEKF for each analysis. Here, the sparse observational
network is used, l2 is ﬁxed to 0. The local patch parameter l is chosen as an optimal settingthat 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 ofserial 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 ﬁxed to4%, the parameter l2 in LEKF is ﬁxed to 0.

3.5.5

Comparison

For comparison of each method, the results of each method are combined and shownin 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 thebest 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 signiﬁcantly
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 thatthe 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 diﬀerent formulation for the polar patches in the LEKF.
As for computational time, Table 3.11 shows timing on a Linux PC with a 2.53GHzIntel 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, especiallyin 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 ﬁeld 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 ﬁeld in the case of the sparse observational
network using 3DVAR (short dashed line), serial EnSRF (solid line) and LEKF (dashedline). The ensemble size for serial EnSRF and LEKF is 30. Parameters are optimal, thatis, 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 withtens 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. Moresatellite 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 observations (26400 observations) are assimilated using a parallel computer with 100 computational 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 α, respectively. The computation is accelerated as
T0 = T ×

1
Nα(3.12)

where T and T 0 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:
58(3.13)

0
TEnSRF= TEnSRF × 5(3.14)0
TLEKF= TLEKF ×

Table 3.12 shows the estimated timimg. By the eﬃcient parallelization and large numberof 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 diﬀerent observational network and investigate the impact of moisture observations. Then, we investigate the impact of vertical error correlations and a diﬀerent
type of error covariance localization. Finally, we apply random perturbation addition toensemble members. The sensitivity results also suggest possible ways to improve the ﬁlter
performance.

1003.6.1

Response with a diﬀerent observational network

So far, we applied mostly the sparse observational network, and we saw lattice-like patterns
in the analysis error ﬁelds. Since the observations were available uniformly in space, theerror ﬁelds 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
ﬁelds. 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 diﬀerent
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 ﬁeld in the cases of 3DVAR (shortdashed 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 ﬁelds
of 3DVAR (top panel) and serial EnSRF (bottom panel). Both ﬁgures have a similar
structure with diﬀerent 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 Paciﬁc 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 eﬀect 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 ﬁelds look more realistic; data-rich regions
101

Figure 3.29: Analysis RMSE of 500hPa height ﬁeld 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ﬁelds in the cases of 3DVAR and serial EnSRF, respectively. In all ﬁgures, it is clearthat the southern hemisphere has much larger errors than the northern hemisphere. By
comparing 3DVAR and serial EnSRF, both ﬁgures show very similar structure with muchlarger amplitudes in 3DVAR. However, there is a structural diﬀerence 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 diﬀerence in the error amplitudes. This is clearer
in the temperature ﬁeld (bottom panel). Although many observations are available in
the northern hemisphere, there are data-poor regions over the Atlantic and Paciﬁc ocean.
Since 3DVAR is more sensitive to the lack of data, 3DVAR tends to be aﬀected more overthe data-poor regions. This could explain why 3DVAR shows larger errors in the northernhemisphere. RMSE of temperature ﬁelds 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 ﬁelds are strongly related to highly nonlinear and complicated physical processes,
moisture ﬁelds tend to show a noisy structure compared to other variables. Since themoisture ﬁelds represent a complicated and noisy physical processes, it is not clear if the
moisture observations actually provide useful information in data assimilation, especiallyin 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 ﬁelds (in meters) inthe 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 is3.0. Dots indicate the observational locations.

104

Figure 3.31: Zonal structure of analysis RMSE of u-wind ﬁeld (top panel) and temperatureﬁeld (bottom panel) in the case of 3DVAR with the realistic observational network. Shades
and contours show analysis RMSE and mean analysis ﬁelds, respectively.

105

Figure 3.32: The same ﬁgure 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 conﬁrm ourresults are not strongly aﬀected by the noisy moisture observations.
In order to investigate the impact of moisture observations, we perform data assimilation cycle experiments excluding moisture observations and compare with the previous
results including moisture observations. Fig.3.33 shows analysis RMSE using serial EnSRF with 10 ensemble members with and without moisture observations. For the moisturevariable (speciﬁc humidity, q), we obtain large improvement by moisture observations as
expected. Importantly, we see positive impact by the moisture observations for othervariables. We also see improvements by moisture observations using LEKF, supporting
moisture observations have positive impact. Thus, it is reasonable to include moistureobservations in our experiments, suggesting our results are not adversely aﬀected by the
moisture observations.
Note that the SPEEDY model has only very simple physics, whereas the physicalprocesses in more sophisticated models are much more complicated. As a result, it is
not straightforward to determine in the operational context if moisture observations causesuch 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 important 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 ﬁeld (a), hight ﬁeld (b), temperature ﬁeld (c)and speciﬁc humidity ﬁeld (d) in the cases of serial EnSRF with 10 ensemble memberswith 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 levelis 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 ensemble members. The local patch parameters are ﬁxed 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 correlation indeed shows a positive impact. With 20 ensemble members, however, it has slightnegative 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 diﬀerence is not statistically signiﬁcant.
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 threedimensional 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 fullvertical 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 is109

Figure 3.34: Impact of the vertical correlation in the analysis RMSE of 500hPa height ﬁeld
using LEKF with 10 (black lines) and 20 (blue lines) ensemble members. The black solidline 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 andwithout the vertical correlation, respectively.

110

Figure 3.35: The same ﬁgure 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 thevertical correlation provides useful information, but sampling errors are also introduced,
even with 20 ensemble members. Thus, the positive and negative eﬀects are cancelled out
and the net eﬀect 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 conﬁguration. However, LEKF
with 10 ensemble members has performed worse than 3DVAR. The diﬀerence between
EnSRF and LEKF is in the way of localization. Serial EnSRF localizes error covariancearound observations with Gaussian-like weighting, LEKF localizes around an analyzed
point with Heaviside-like step weighting. Thus, we applied observational error covariance111Localization scale (σ)

∞

1.0

√
1/ 2

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 ﬁeld using LEKF with 10 ensemble memberswith 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 ﬁeld using LEKF with 10 ensemble members and local patch parameter l = 1 with various localization scales of observational error covariance localization. We can see a large improvement of LEKF withthe observational error covariance localization in the case of 10 ensemble members. If wechoose 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 whenlocal 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 thelocalization when the ﬁlter 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 ﬁelds with
diﬀerent localization settings. After the initial one-month transient period, the three lines112

Figure 3.36: Analysis RMSE of 500hPa height ﬁeld using LEKF with 10 ensemble members. Here, observational error covariance localization is applied. The solid, long-dashedand 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 thevery initial transient period (ﬁrst 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 eﬀect. 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 modiﬁed 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 diﬀerent localization
weighting function turned out to be very small. When we took 1 × 1 squares aroundobservations, 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 scaleof 1. Here, we used 10 ensemble members.
In EnKF formulations, we applied globally the same localization length scale ingrid 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 perturbations 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 foundthat the augmentation of the background error covariance with BVs reduced analysis errors 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) investigated the eﬀect of the ”stochastic seeding” in a breeding cycle on the Lorenz-96 model,

114

Figure 3.37: Analysis RMSE of 500hPa height ﬁeld using LEKF with 10 ensemble members. 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 ignoredin 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 bepossible that ensemble members lose signals of instabilities.
Expecting that the random perturbation addition or ”stochastic seeding” has similareﬀects 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 ﬁeld using LEKF with and without random perturbation addition. With
random perturbation addition, the RMSE is smaller almost everywhere or at least is notworse than the case without random perturbation addition.

115

3.6.6

Summary

In this section, we conﬁrmed that our results in the previous section were reasonably
robust with the experimental settings. Using an irregularly distributed realistic observational network, both 3DVAR and serial EnSRF show larger errors than the regularly
distributed sparse network. 3DVAR is more strongly aﬀected by the irregular distributionof 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, accounting for local vertical correlation improve the results, but full vertical correlation worsens
the ﬁlter. With 20 ensemble members, the eﬀect of considering vertical correlations is
not signiﬁcant. These results suggest that the vertical error independence is a reasonablechoice in the SPEEDY model. LEKF shows signiﬁcant improvements by the observational 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 nothave 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 analyze the structure of analysis errors (which are not known in real data assimilation) andto 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 aquasi-geostrophic model (Kalnay et al. 2005). We compare forecast errors and analysis
increments for 3DVAR and EnKF.
Fig.3.38 shows the background error ﬁelds (forecast minus truth, shaded) and theanalysis increment ﬁelds (forecast minus analysis, contour) at an arbitrary time in the caseof 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 ﬁeldof EnKF captures ﬂow-dependent structures of the background error ﬁeld, 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 likecircles 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 ﬁeld 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 corresponding to observational location. Although correlations between the analysis increment andbackground error are small (around 0.3) both in 3DVAR and EnKF, EnKF shows a clear
advantage because of better background ﬁelds.

117

Figure 3.38: Background error ﬁeld (shaded) and analysis increment ﬁeld (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 ﬁelds 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.2Analysis 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 scaleof 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 tomeasure the perturbation amplitude:
Ã

1
cp
E=
u2 + v 2 + T 2 + Rd Tr
2
Tr

µ

ps
p0

¶2 !

(3.15)

where E is energy, and constant-pressure heat capacity cp = 1005.7Jkg−1 K−1 , dry-air gas
constatnt Rd = 287Jkg−1 K−1 , reference temperature Tr = 270K, and reference pressurep0 = 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 errors. ”SPREAD” shows forecast/analysis ensemble perturbation spreads, indicating the
growth of ensemble perturbations. The perturbation amplitudes are averaged in timeusing whole February (112 samples) so that the ﬁgure 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 ﬁrst few
days. However, in their EnKF scheme, they added a model error covariance term to each120

SPEEDY-3DVAR/EnSRF 5-DAY FORECAST
40
RMSE(3DVAR)SPREAD(EnSRF)
RMSE(EnSRF)

35

2 2

ENERGY [m /s ]

30

25

20

15

10

5

0
0

24

48

72

96

120

TIME [HOUR]

Figure 3.40:

The growth of the initial analysis errors and ensemble perturbationsin 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 thegrowth 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 thestatistical average and are unrelated to the daily ﬂow. Thus, we expect these perturbations 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 arealmost parallel. Thus, ensemble perturbation and analysis error growths show very similar
growth rates (doubling in about three days), which implies both perturbation ﬁelds 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 conﬁguration, it seems 4% inﬂation is too large,
which is why the ensemble spread is larger. To conﬁrm 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, theanalysis 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 ﬁelds and ensemble spread

Since EnKF works as expected, ensemble members are expected to capture the error
structures. Thus, we plot the error ﬁelds and ensemble spread in the case of serial EnSRF
with 30 ensemble members and the Schur product length scale of 3.0. The top panelof Fig.3.41 shows the analysis error ﬁeld at an arbitrary time, and the bottom panel
shows the forecast error ﬁeld at the same time. The analysis and forecast error ﬁelds aredeﬁned as analysis/forecast minus the true ﬁeld. The shaded ﬁeld in Fig.3.41 shows theanalysis/forecast error ﬁeld, 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 theerror ﬁeld and the ensemble spread. However, the correspondence is not very good, thepattern correlations between error ﬁelds and ensemble spread are 0.25 in the analysis and0.18 in the 6-hour forecast. By contrast, the pattern correlation between analysis andforecast ensemble spreads are very high (0.96). Similarly, the pattern correlation between
analysis and forecast error ﬁelds are also high (0.74).

3.7.4

Forecast errors, analysis errors, and bred vectors

In order to see the shape of the error ﬁelds, we plot the 2-day forecast error and analysis
error ﬁelds valid at the same time. Fig.3.42 shows 2-day forecast error ﬁeld (shaded) andanalysis error ﬁeld (contour) of u-wind at the 4th level (σ = 0.51) at an arbitrary time inthe case of serial EnSRF with 10 ensemble members and the localization scale of 1. It isclear that the two ﬁelds have extremely similar shapes, conﬁrmed by their high patterncorrelation of 0.76. 2-day forecast error ﬁeld is regarded as a dynamically growing error,
the analysis ﬁeld contains growing errors. Even with EnKF that considers dynamical error

123

Figure 3.41: The error ﬁelds (ananlysis/forecast minus truth; shaded) and the ensemble
spread of 500hPa height ﬁeld (contour) at an arbitrary time. The top and bottom panels
show analysis and background (6-hour forecast), respectively, both at the same time. Someparts have similar structures indicated by rectangles.

124Figure 3.42: 2-day forecast error ﬁeld (shaded) and analysis error ﬁeld (contour) of u-windat 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
ﬁelds are 0.76.

structure, the analysis ﬁelds still contain growing errors. Fig.3.43 shows the same ﬁgurein the case of 30 ensemble members, which still shows high correlation (0.79). Fig.3.44
shows the same ﬁgure in the case of 3DVAR, which also shows high correlation (0.74).
The three ﬁgures show the same time with the same truth, so it is expected the
dynamically growing error has similar shape. However, all the three ﬁgures show diﬀerent
shapes, which implies there are many possible dynamical ”errors of the day” rather thana dominance by bred vectors.
One can observe high correlation of the two ﬁelds not only at this moment for the
speciﬁc variable, but also at every time for all variables. Fig.3.45 shows temporal meanof the correlation coeﬃcients for four model variables (u, v, T , q), all of which show high125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 ﬁelds are 0.79.

Figure 3.44: The same as Fig.3.42, but in the case of 3DVAR. The correlation between
the two ﬁelds are 0.74.

126Figure 3.45: Temporal mean of correlation coeﬃcients between 2-day forecast error andanalysis error ﬁelds 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 ﬁelds, 5-day and 2-day forecast
errors, analysis error, E-dimension, and bred vector (BV) ﬁelds are shown in Fig.3.46.
E-dimension, initially proposed by Patil et al. (2001), is a measure of substantial dimensionality of the space spanned by ensemble members, the deﬁnition 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 ﬁeld consisting of random numbers. The time intervalof 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-variablecorrelations, (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 ofPaciﬁc ocean, western Siberia, and the southern end of south America. We ﬁnd 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 ﬁelds (e) and (f), we seea similar signal in the western Siberia with smaller amplitudes. The growing errors shown
in BV and forecast ﬁelds almost disappear in analysis error ﬁelds.

3.7.5

E-dimension and explained variance

Szunyogh et al. (2005) have shown a strong anti-correlation between E-dimension andexplained variance. E-dimension φ is a measure of substantial dimensionality of the spacespanned by ensemble members, deﬁned asP( m λi )2
φ = Pi=1m
2
i=1 λi(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 deﬁned as

ξ=||GG> (xf − xt )||2||(xf − xt )||2

(3.17)The top and bottom panels of ﬁg.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, Edimension is computed using a 3 × 3 local patch for all variables, inter-variable correlations 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 4thlevel 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 etal. (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 thescatter plot, indicating a high anti-correlation of -0.66. The zonal mean is not essential toobtain the high anti-correlation between E-dimension and explained variance.
However, if we ﬁx the time, we do not observe the strong anti-correlation. Fig.3.51
shows the same ﬁgure as Fig.3.49 but at a ﬁxed time. Fig.3.52 shows the scatter plot between E-dimension and explained variance shown in Fig.3.51, indicating low anti-correlationof -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 SouthAmerica, and northern Paciﬁc, we ﬁnd correspondence of low E-dimension and large explained variance. In these several parts, EnKF could capture the true error shape andcontribute 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 ﬁeld of EnKF captures ﬂow-dependent structures of the
background error ﬁeld, 3DVAR cannot capture the ﬂow-dependent structures at all.
• 3DVAR showed a striped pattern in error ﬁelds, 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 thegrowth rate of the analysis error.
• Error and ensemble spread ﬁelds show similar signals in a few areas, but with low
pattern correlations of about 0.2.
• 2-day forecast error and analysis error ﬁelds have very similar shapes, with correlations about 0.75. However, the error ﬁelds are diﬀerent with diﬀerent data
assimilation schemes, suggesting the possibility of multiple ”errors of the day”.
• 5-day forecast error, 2-day forecast error and bred vector ﬁelds show similar signals
in low E-dimensional areas.
• Time-mean ﬁelds of explained variance and E-dimension are anti-correlated with pattern 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 forone month (112 samples).

133

Figure 3.50: Scatter plot of the time mean E-dimension and explained variance at the 4thlevel, 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 andapplied on the SPEEDY model under the perfect model assumption. Our results showthat the EnKF clearly outperforms 3DVAR. One surprising result is that for diﬀerent
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 diﬀerent
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 ﬁelds 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 anarbitrary time, indicating anti-correlation -0.26.

Overall, our results suggest that serial EnSRF outperforms LEKF, but this advantage 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 onl2 . A natural way is to use Gaussian-like weighting according to the distance from thecenter 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 overthe 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 lesscomputational 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) areapplied on the SPEEDY model. EnKF clearly outperforms 3DVAR under the perfectmodel assumption. In this chapter, we remove the perfect model assumption using the
NCEP/NCAR reanalysis (NNR) ﬁelds (Kalnay 1996) as the nature run. We explore thefollowing 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 eﬀects 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 diﬀerences between
NNR ﬁelds and the SPEEDY model 6-hour forecasts initialized by NNR ﬁelds. Then, theeﬀects 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 ﬁelds. We describe theory on model bias estimation following Deeand da Silva (1998) and apply the estimation method with EnKF. We also apply low-orderbias estimation, assuming the bias ﬁeld is spanned by a limited number of base ﬁelds. Asbase ﬁelds, we use the model error bias ﬁeld and EOFs.
138

4.1.2OSSE 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 forboth generating observational data and performing data assimilation. This perfect modelexperiment is over optimistic partly because such observational data do not contain theerror of representativeness. Alternatively, Atlas (1997) suggested to generate observations
using a more sophisticated and higher resolution model, diﬀerent from the model used indata assimilation. This considers the error of representativeness as well as the imperfectionof the model. Similarly in this chapter, we use NNR ﬁelds to generate observational data,
the NNR ﬁelds are generated by a relatively recent sophisticated model (1995 NCEP
model) and data assimilation system. The diﬀerence from the Atlas’s approach is that
we consider the real nature. Real observations are assimilated in NNR ﬁelds, NNR ﬁelds
simulate the evolution by the unknown ”true nature” model.

4.1.3

NCEP/NCAR reanalysis (NNR)

The NCEP/NCAR reanalysis (NNR) ﬁelds (Kalnay 1996; Kistler et al. 2001) are a fairly
good estimate of the true atmospheric state for the past 50 years. They are generated bya relatively recent (1995) sophisticated model and data assimilation system. We assume
the NNR ﬁelds as a proxy for the truth. The validity of this assumption is beyond thescope of this research. We randomly perturb the NNR ﬁelds to simulate observations withthe same standard deviation (Table 3.2), we verify our analysis ﬁelds against the NNR
ﬁelds.
NNR data are freely available on the Internet. The NNR data on pressure levels

139have a diﬀerent grid structure (144 × 73 × 17) than the SPEEDY model (96 × 48 × 7). The
NNR ﬁelds need to be interpolated to obtain the ﬁelds on the SPEEDY model grids. The
horizontal interpolation is a ﬁrst order linear interpolation. No vertical interpolation isrequired to obtain pressure level data because the reanalysis data have all seven pressurelevels: 925, 850, 700, 500, 300, 200, and 100hPa. A ﬁrst order linear interpolation isapplied 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 ﬁeld is interpolated for the target level. log p is used as a vertical coordinate
for the linear interpolation.
Since the original NNR ﬁelds have a higher resolution, they contain wavenumber
components higher than those which the SPEEDY model can resolve. Wavenumber components higher than 30 cannot be detected in the SPEEDY model since it has spectral
triangular truncation of 30 wavenumbers. The reanalysis ﬁelds are smoothed out by thespectral transformation of the SPEEDY model.
For the present experiments, we downloaded the NNR ﬁelds from the NCAR’s website and generated the interpolated and smoothed ﬁelds 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 ﬁelds ”reanalysis ﬁelds”.

4.2 Model error statistics
In this section, model error statistics (bias and EOFs) are computed and shown.

140

4.2.1Model error biases

In order to estimate the model error against the reanalysis ﬁelds, the reanalysis ﬁelds are
given as initial conditions to the SPEEDY model. Since the reanalysis ﬁelds are evolved
by an unknown model expected to simulate the true atmospheric dynamics, the diﬀerences
between the SPEEDY model 6-hour forecasts and reanalysis ﬁelds valid at the same time
represent the model errors. We assume the diﬀerence is a sample of model errors.
Fig.4.1 shows the diﬀerences of the three prognostic variables (u, v and T ) between
the SPEEDY model 6-hour forecasts and reanalysis ﬁelds, temporally averaged for twomonths using 235 samples in the period from January 1, 1982 to February 28, 1982.
Temporally averaged diﬀerences are interpreted as the model bias. Large areas are notshaded by colors, suggesting that the SPEEDY model simulates the real atmosphere fairly
well. The largest model bias of the wind ﬁelds can be seen in polar regions. Somediﬀerences are caused by the orographic eﬀect, most clearly seen in the u-wind ﬁeld at
200hPa level, where there is a clear negative bias downstream of the Himalayas. As
for temperature ﬁelds, 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 ﬁgure as Fig.4.1 but for the variables q (speciﬁc 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 aﬀected by the psbias ﬁeld. It seems that the orographic diﬀerence in the reanalysis ﬁeld and the SPEEDY
model aﬀects the ps model bias. However, it still shows some biases over the ocean, e.g.
positive bias over the western tropical paciﬁc basin, which disappears in the z bias ﬁeld
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 ﬁelds.

142

Figure 4.2: SPEEDY model biases of q, z [m] and ps [Pa] relative to the NCEP/NCAR
reanalysis ﬁelds. Each vertical level has diﬀerent scales for q ﬁelds, 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 deﬁned as eigenvectors of the covariancematrix 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 deﬁned by subtracting the time-mean model bias ﬁeldx̄(s) from time series of model error ﬁelds 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 composedof eigenvalues, where N denotes the number of the space index. Usually, N is a numberof O(104 ) or larger, the number of time index m is smaller. As a result, the covariancematrix 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/2 V>(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 satisﬁed. Using eq.(4.3), we
get EOFs as follows:144x0 (t, s)> x0 (t, s) =³

UD1/2 V>´> ³

UD1/2 V>´

= VD1/2 U> UD1/2 V>= VDV>(4.4)

V is equal to that in eq.(4.1), containing EOFs. There is an algorithm to ﬁnd the largest
singular values and corresponding singular vectors, we obtain ﬁrst few EOFs with a little
computation.
Here, we take space index s by combining variables (u, v, T , and q) but separatinglevels. In this way, inter-variable correlations are considered but no correlations between
levels are considered. Fig.4.3 shows the ﬁrst EOF (top panel) and the second EOF (bottom panel) of temperature ﬁelds at the bottom level (shades) and surface pressure ﬁelds(contour). There is a clear signal corresponding to diurnal variation with wave number 1
over longitudes. We also see the phase diﬀerence between the ﬁrst and second EOFs. The
SPEEDY model considers only daily updated forcings and no diurnal variation. Temperature at the bottom level is strongly aﬀected 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 ﬁrst four EOFs. Fig.4.4 shows timeseries of expansion coeﬃcients of the ﬁrst four EOFs for a week. We see a clear signal ofthe diurnal cycle for the ﬁrst and second EOFs. The second EOF shows diurnal ﬂuctuation with a slightly smaller amplitude and a diﬀerent phase than the ﬁrst EOF. Fig.4.5
shows time series of expansion coeﬃcients of the ﬁrst EOF (top panel) and the third andfourth EOFs (bottom panel) in a longer time range. Diurnal ﬂuctuation is dominant inthe ﬁrst 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 ina 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ﬁelds themselves and expansion coeﬃcients indicate that with this short training period,
only the ﬁrst two EOFs (representing the error due to the absence of diurnal cycle) canbe estimated in a statistically signiﬁcant fashion.

4.3 Eﬀects of model errors on data assimilation
While the previous section showed bias and EOFs of the model errors, the eﬀects of systematic model errors on data assimilation are not clear. To see the eﬀects, we perform thesame data assimilation experiments as the perfect model experiments using the reanalysis
ﬁelds as truth. In other words, the same data assimilation systems are applied, the only
diﬀerence is the nature run is replaced by the reanalysis ﬁelds. 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 ﬁeld. Even EnKF with 30 ensemble members, the best among the three, shows about 60 meters of analysis RMSE.
3DVAR shows the largest temporal ﬂuctuation and largest mean errors. However, theperformance of EnKF is more than 10 times worse than in the perfect model case, and theratio of 3DVAR to EnKF is clearly smaller than in the perfect model case. This meansthat EnKF is more strongly aﬀected 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 ﬁelds at the bottom level (shades) andsurface pressure ﬁelds (contour) of the ﬁrst and second EOFs, respectively.

147

150
EOF1
EOF2
EOF3
EOF4

100

EXPANSION COEFFICIENTS

50

0

−50

−100

−150

−200

1

2

3

4

5

6

7

8

TIME[DAY]

Figure 4.4: Time series of expansion coeﬃcients of the ﬁrst to fourth EOFs for a week.
Solid, dashed, dotted, dash-dotted lines show ﬁrst, second, third, and fourth EOFs, respectively.

148

150
EOF1

100

EXPANSION COEFFICIENTS

50

0

−50

−100

−150

−200

5

10

15

20

25
30
TIME[DAYS]

35

40

45

50

150
EOF3
EOF4
100

EXPANSION COEFFICIENTS

50

0

−50

−100

−150

−200

5

10

15

20

25
30
TIME[DAYS]

35

40

45

50

Figure 4.5: Time series of expansion coeﬃcients of the ﬁrst 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 ﬁeld of 3DVAR (short dashed line) and
EnKF with 10 ensemble members (solid line) and 30 ensemble members (long dashedline) in the case of the regular sparse observational network (Fig.3.10) in the presence ofmodel errors.

150

Figure 4.7: Analysis RMSE of 500hPa height ﬁeld of 3DVAR and EnKF with 30 ensemble
members in the case of the sparse and realistic observational networks in the presence ofmodel 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 dottedline 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 ﬁeld in the cases of both the sparse and realistic observational networks. Clearly,
the diﬀerence between 3DVAR and EnKF is larger in the case of the realistic observational
network. 3DVAR is strongly aﬀected by the diﬀerent 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. Signiﬁcant diﬀerenceis only seen in height and temperature ﬁelds at a few levels; EnKF shows advantage inheight ﬁelds at middle and upper levels, 3DVAR shows advantage in temperature ﬁelds
at 500, 200, and 100hPa levels. The fact that EnKF is less sensitive to the choice ofobservational network than 3DVAR is generally valid for other variables at other levels
except temperature and speciﬁc humidity, and height ﬁelds at upper levels. There isalmost no diﬀerence among the ﬁve data assimilation methods in speciﬁc humidity ﬁelds.

4.4 Model bias estimation
The previous section showed that model errors aﬀect data assimilation, especially in EnKF.
The advantage of EnKF to 3DVAR became smaller than the perfect model case. In thissection, model bias estimation (Dee and da Silva 1998) is introduced and applied to
improve the ﬁlter performance in the presence of model errors.

4.4.1

TheoryWe 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 amodel bias:
x̃fi+1 = M (x̃ai ) − bai(4.5)

Here, x̃ denotes an unbiased state. Since the model M has a bias b, it is necessary tosubtract the bias to get the unbiased forecast. The bias also evolves in time:bfi = Mb (bai−1 )

152

(4.6)

Figure 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-dashedlines), 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 realistic observational network, respectively. The four panels (a),(b), (c) and (d) correspond to u-wind, height, temperature (T ) and speciﬁc humidity (q)ﬁelds, respectively. The observational error standard deviations are shown as thin solidlines 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 ﬁltering equations as in Section
2.2, thus,x̃ai = x̃fi + K(yo − H x̃fi )(4.7)

where the optimal weight K is given as
K = Pf H> (HPf H> + R)−1

(4.8)

cf. eq.(2.15). Here, the error covariance matrix Pf is deﬁned for the unbiased variable asD

Pf = (x̃f − xt )(x̃f − 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 theensemble forecast, and the ensemble update is given by the SRF (square root ﬁltering) tosatisfy eq.(4.10).
Similarly for the bias, we deﬁne 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 deﬁne the diﬀerences from the true bias bt :
δba = ba − bt

(4.12)= bf − bt(4.13)

δbf154

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δx̃f + Hδbf´

(4.14)

where δxf is deﬁned as the same as eq.(2.3) in Section 2.2, and we assumed unbiased
forecast error:
δx̃f= xf − xt − bt

(4.15)= x̃f − xt

(4.16)= δxf − bt

(4.17)

Bias analysis error covariance Pab is written
Pab =

D

δba (δba )>

¿³

=
¿³

=

E
´³

δb + Kb (δy − Hδx̃ + Hδb )
f

o

f

f

(I − Kb H)δbf − Kb δyo + Kb Hδx̃f´> Àδb + Kb (δy − Hδx̃ + Hδb )
f

´³

o

f

f

(I − Kb H)δbf − Kb δyo + Kb Hδx̃ff > >= (I − Kb H)Pfb (I − Kb H)> + Kb RK>
b + Kb HP H Kb´> À

(4.18)

where bias forecast error covariance is deﬁned asD

Pfb = δbf (δbf )>E(4.19)

and all the cross covariances are assumed to be zero:Dδx̃f (δyo )>Dδbf (δyo )>Dδx̃f (δbf )>
155

E

= 0

(4.20)= 0(4.21)= 0(4.22)EE

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,∂(trace(Pab )) = 0∂Kb(4.23)

is satisﬁed (cf. eq.(2.11)). Using the formula (2.12) and (2.13), we get the optimal weight
Kb = Pfb H> (HPfb H> + HPf H> + 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 − Kb H)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) isjust
bfi = bai−1(4.27)In this case, the bias has no dynamics, and the low dimensional assumption may not bereasonable.
156

Alternatively, we can choose a limited number, say k, of base ﬁelds to span the biasb (Dee and da Silva 1998; Danforth, Kalnay, pers. comm.):b = Tβ(4.28)

Here, T is an N × k matrix where each column denotes base ﬁelds, and β denotes a kdimensional vector. In this way, β is the variable to be estimated, T is ﬁxed. For the best
estimate of β, we substitute eq.(4.28) into the previous equations derived for b. Here, theerror covariance of b is substituted asDδb(δb)>

Pb =

E

D

=

Tδβ(Tδβ)>DEE= T δβ(δβ)> T>= TPβ T>

(4.29)

Thus, substituting eqs.(4.28) and (4.29), we construct the Kalman ﬁltering equations forβ as follows:βifa= βi−1(4.30)

Pfβ = (1 + δ)Paβ(4.31)

βia = βif − Kβ (yo − Hxfi + HTβ fi )(4.32)Kβ = Pfβ T> H> (HTPfβ T> H> + HPf H> + R)−1

(4.33)Paβ = (I − Kβ HT)Pfβ(4.34)

where δ denotes a covariance inﬂation parameter for Pβ . If the dimension of β is smallenough, 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 formulationof Kalman ﬁltering. We use the serial EnSRF system with minimal modiﬁcations. The
diﬀerence is in the gain computation. Since we need to consider not only Pβ 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 modiﬁcation is not complicated.
The initial bias ensemble is generated using Gaussian random numbers. A random
perturbation ﬁeld is constructed by putting a Gaussian random number at each grid point
and smoothing high frequency components using a Lanczos ﬁlter. The amplitudes arenormalized using observational error size, so that each random variable has an appropriate
order of magnitude. We used 10 ensemble members, where 5 random ﬁelds 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 ﬁlter was
unstable and diverged rapidly. In order for the ﬁlter to be more stable in the initialtransient period, the bias correction eq.(4.5) is modiﬁed asx̃fi+1 = M (x̃ai ) − αbai(4.35)

where α is a scalar parameter that weights the bias correction. Even with α = 0.01, one
percent bias correction, we still observed unstable ﬁltering. Fig.4.9 shows the analysis
158

Figure 4.9: Analysis RMSE of 500hPa height ﬁeld of EnKF with 10 ensemble members in
the case of the sparse observational networks in the presence of model errors. The dashedand solid lines show the cases with and without Dee and da Silva’s full dimensional bias
estimation, respectively.

RMSE of 500hPa height ﬁeld with/without bias estimation. Here, we used α = 0.01. The
time period ﬁnishes before February, because the ﬁlter diverges just after the ﬁnal time of
the ﬁgure even with α = 0.01. In the ﬁrst 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 ﬁeld
cannot have dynamically growing modes. Therefore, the low-dimensional assumption for159

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 ofobservations, larger number of observations helps the ﬁlter converge faster. With fewer
observations and without dynamics, the initial spin-up period of the ﬁltering becomes
longer.
Since the full dimensional estimation made the ﬁlter unstable, we applied lowdimensional bias variable estimation. We used the mean bias ﬁelds 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 β
is ﬁxed 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 ﬁeld. Thin red line shows
6-hour forecast errors initiated from the reanalysis ﬁelds. These forecasts provide background errors purely caused by model errors. Clearly, EnKF shows smaller analysis RMSE
with bias estimation. There is no signiﬁcant diﬀerence between low order bias estimation(green short-dashed line) and constant bias subtraction (blue dotted line). 3DVAR (bluelines) shows slightly better performance with constant bias subtraction, but the improvement is smaller than with EnKF. Fig.4.11 shows the values of all the seven componentsof 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 theconstant bias subtraction shows similar performance.
Fig.4.12 shows the spatial distribution of analysis RMSE of 500hPa height ﬁelds
with and without bias estimation. Clearly, the amplitude is reduced more in the case withbias estimation (bottom panel) than the case without bias estimation (top panel). The
large diﬀerence can be seen for example over the Andes of South America, Himalayas,

160

Figure 4.10: Analysis RMSE of 500hPa height ﬁeld for 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
solid and long-dashed lines of Fig.4.6. Thin red line shows 6-hour forecast errors initiatedfrom the reanalysis ﬁelds.

161

SPEEDY-EnSRF BIAS VARIABLES
1.5

1

0.5

1
2
3
4
5
6
7

0

-0.5
40

80

120

160

200

TIME [DA CYCLES]

Figure 4.11: Values of all the seven components of the bias variable when the mean
bias ﬁeld 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 thebottom 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 ﬁelds are shown in
Fig.4.13. The case without bias estimation (top panel) shows much larger amplitudes thanthe case with bias estimation (bottom panel), large error spots exist in similar locations.
If we subtract the error bias from the error ﬁelds and take the RMS, i.e. analysis error
standard deviation, the diﬀerence between the cases with and without bias estimation
becomes small. Fig.4.14 shows that the diﬀerence between the top and bottom panels are
smaller than the previous ﬁgures. 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 difference is 30-40 meters in height ﬁelds, 2-3 m/s in wind ﬁelds, 1-3 K in temperatureﬁelds, all of which are fairly large improvements. There is almost no diﬀerence between
the cases with low order bias estimation (green short-dashed line) and with constant bias
subtraction (black dotted line), only speciﬁc humidity RMSE shows signiﬁcant diﬀerence.
3DVAR and EnKF show similar performance without bias subtraction (solid lines). However, 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 speciﬁc humidity, the situation is opposite, 3DVAR shows smaller RMSE than EnKF.
With low order bias estimation, EnKF shows almost identical performance as 3DVARwith constant bias subtraction. Concentrating on dynamics, we obtain much better per-163

Figure 4.12: Spatial distribution of analysis RMSE of 500hPa height ﬁelds (shades in meters) in the case of imperfect model experiments using EnKF with 10 ensemble members.
Contours show mean analysis ﬁelds (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 ﬁelds of 500hPa height (shades in meters) in the caseof imperfect model experiments using EnKF with 10 ensemble members. Contours show
mean analysis ﬁelds (in meters). Top and bottom panels show the cases without and withbias estimation, respectively. Dots indicate the observational locations.

165

Figure 4.14: Analysis error standard deviation ﬁelds of 500hPa height (in meters) in thecase of imperfect model experiments using EnKF with 10 ensemble members. Contours
show mean analysis ﬁelds (in meters). Top and bottom panels show the cases without andwith bias estimation, respectively. Dots indicate the observational locations.

166

formance using EnKF with bias estimation and constant bias subtraction in wind andheight ﬁelds. However, for variables strongly related to physical processes, the advantage
of EnKF to 3DVAR is less signiﬁcant. There is signiﬁcant improvement of humidity ﬁelds
with low order bias estimation compared to constant bias subtraction, similar to 3DVARwith constant bias subtraction.
Fig.4.16 shows the zonal structure of analysis RMSE of u-wind and temperatureﬁelds in the case without bias estimation. Fig.4.17 shows the same ﬁgure but in the casewith 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 ﬁelds show aclear striped pattern according to the observational locations. We see a large reductionof temperature errors in the lower altitudes in the northern hemisphere in the case withbias estimation.
We include the ﬁrst few EOFs as bases to span the bias variables. Table 4.1 shows
analysis RMSE of 500 hPa height ﬁeld with various numbers of EOFs. The mean bias
ﬁeld and ﬁrst two EOFs give the best performance, in agreement with the argument in
Section 4.2.2 the higher order EOFs do not contain statistically signiﬁcant information
due to the shortness of the training period. Since the diurnal cycle components and thebias are well deﬁned, even with a short training period, their online correction leads tosigniﬁcant 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 ﬁelds, as suggested

167

Figure 4.15: Analysis RMSE at the all pressure levels temporally averaged for one monthafter the initial one-month spin-up period (112 samples) in the cases of EnKF with loworder 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 ) andspeciﬁc humidity (q) ﬁelds, respectively. The observational error standard deviations areshown as thin solid lines wherever applicable.

168

Figure 4.16: Zonal structure of analysis RMSE of u-wind ﬁeld (top panel in m/s) andtemperature ﬁeld (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 ﬁelds, respectively.

169

Figure 4.17: The same ﬁgure 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 ﬁeld using EnKF with 10 ensemble members
and low order bias estimation using various numbers of bias basic ﬁelds. ”Uncorrected”
denotes the case without bias correction. ”Mean bias”, ”EOF1”, and ”EOF2” denotes thecases with low-order bias correction using the mean bias ﬁeld, using the mean bias ﬁeld
and ﬁrst EOF, and using the mean bias ﬁeld and ﬁrst 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 diﬀerence between thereanalysis ﬁelds and 6-hour forecasts started from the reanalysis ﬁelds. We computed themean bias and EOFs from the model error samples, and used them as bases of modelbias. Then, we found the low order bias estimation has a positive impact even with asingle statistical bias ﬁeld. Constant bias subtraction has also shown similar performance,
although it resulted in larger RMSE in speciﬁc humidity ﬁelds. 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 speciﬁc humidity. When we applied low order bias estimation within EnKF,
EnKF outperformed 3DVAR everywhere and for all ﬁelds. It is important to note that
we used the dependent sample for the statistics, suggesting we may have overestimatedthe performance. It is expected that if we sample with other terms or years to obtain thebasic bias ﬁelds, our results would be worse than what we show in the present research.
However, other experiments have shown similar basic bias ﬁelds even using sampling of171

other years (Danforth et al. 2005).
Another problem is that the SPEEDY model has much larger errors than sophisticated operational models. We used the NNR ﬁelds as a best estimate of the true atmosphere and took the diﬀerence from the SPEEDY model forecast to obtain samples ofmodel 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 andforecasts, the diﬀerences between the forecasts and analyses by the same model, i.e. theanalysis increments, will be samples of model errors. The model bias ﬁeld obtained fromthe analysis increments was very diﬀerent from the one obtained from the NNR ﬁelds.
Using the analysis increments as model error samples, we obtained worse results (notshown). 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 ﬁelds such as
NNR and ECMWF (European Centre on Medium-range Weather Forecasting) reanalysis
ﬁelds.

172

Chapter 5
Summary and future directions
5.1 Summary
Our goal was to develop a path towards the operational application of ensemble Kalmanﬁltering (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 eﬀects 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 ﬁlter (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, Ottet 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 theone 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ﬂationthat 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 (Molteni2003). The SPEEDY model is a fast but relatively realistic model allowing a comparison of
methods that address the ﬁrst 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 diﬀerent methods except in regions identiﬁed as having low ensemble dimension (Edim), where they are also similar to bred vectors. Our results suggest that for the SPEEDY
model serial EnSRF outperforms LEKF, but that their diﬀerence is substantially reduced if
we either localize the LEKF error covariance or increase the ensemble size. Using stochastic
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 ﬁnd that the full
dimensional model bias estimation fails. However, if instead we assume that the bias islow 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-dimensionalvariational method (4DVAR)• Application of four-dimensional EnKF (4D-EnKF, Hunt et al. 2004) to treat asynoptic observations• The eﬀect of using a nonlinear observational operator H• The eﬀect 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 Meteorological Agency (JMA), the UK MetOﬃce (UKMO), and the Canadian MeteorologicalCentre (CMC) applied 4DVAR operationally. However, there is no research thus far comparing EnKF with 4DVAR. There are many asynoptic observations in practice, 4DVAR
can treat them intrinsically. Alternatively, Hunt et al. (2004) proposed 4D-EnKF, assuming ensemble forecasts can propagate observational increments in time. When the175

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 canconsider non-diagonal components of the observational error covariance R with less additional computation than other methods. Since more satellite observations with nonlinear
observational operator and correlated errors were available recently, it is important toinvestigate these eﬀects.
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 investigationon 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 tomodel errors that we did not discuss so far: model error covariance. Ensemble members
are expected to represent error structures, but model error covariance could aﬀect themeven 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 byconsidering 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 inerror representation of an ensemble, improving analysis as well. However, we need toconsider the eﬀect of the improved forecast in the bias estimation, since the forecast
model is diﬀerent 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 ensemble 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 forecasting. Ensemble members generated by EnKF are not positive-negative pairs, Wang et al.
(2004) have shown a centered spherical simplex ensemble provides better ensemble forecasts than an ensemble of positive-negative pairs with a limited ensemble size. The ensemble 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 ﬁelds 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 directly connected to peoples lives. It is important that we apply better methods in regional/mesoscale prediction systems. There are not many results on ensemble forecasting177

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 boundary conditions for regional ensemble members. In addition, the scales of weather phenomena expressed by regional/mesoscale models are diﬀerent from global models. Asynoptic observations capturing such mesoscale phenomena are more important and should betreated 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 dynamical 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 thatthey 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 ﬁrst 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

EntitynxINdimension of the systemmINensemble size

hdfx(m)

IN

HEf

r

IN

observational error variancedxf(nx,m)IN

forecast ensemble perturbations Efgain(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 using eqs.(2.45) and (2.48), that is,
Ea = (I − αKH)Ef(A.2)= Ef − αKHEf(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 productusing 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 inputof K for the subroutine enkf_serial is localized. Analysis equation eq.(2.1), that is,x̄a = x̄f + 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 other180

Argument

IN/OUT

EntitynxINdimension of the systemmINensemble size

hdfx(m)

IN

HEf

r

IN

observational error variancedxf(nx,m)IN

ensemble perturbations Efgain(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 perturbations.

Argument

IN/OUT

Entityscale

INlength scale lc

dist

IN

distance dfactorOUT

Schur product factor S(r)Table A.3: Arguments of the subroutine enkf schur that computes the factor of the Schurproduct for localization.

181

Argument

IN/OUT

Entityh(nobsl,ndiml)INobservational 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 + x̄f in a local patchxa(ndiml,m)OUT

analysis ensemble perturbation Ea in a local patch

xabar(ndiml)

OUT

analysis ensemble mean x̄a

Table 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.

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 deﬁned as global
variables of the module, need to be determined in advance. Table A.4 gives arguments ofthe 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 ﬁlter technique (Purser et al. 2003a)
B.1 Introduction
As introduced in Section 3.2.2, recursive ﬁlter (RF) is applied to realize the Gaussian
shape of the horizontal correlation in 3DVAR. RF is a quasi-Gaussian ﬁlter which has animportant advantage in computational cost.
3DVAR requires forward and adjoint transformations. Thus, given Gaussian transformation 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, and LT gives its adjoint. Although this method is currently
used in some operational 3DVAR systems, this matrix application is computationallyexpensive.
Alternatively, RF is much cheaper. RF consists of a set of forward and backward
processes mutually adjoint:

Bi = βAi +Ci = βBi +

n
X
j=1
n
X

αj Bi−j(B.2)

αj Ci+j

(B.3)

j=1

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 ﬁltered output signal. The subscript i denotes the spaceindex 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 ﬁnd coeﬃcients α and β. This appendix
introduces RF theory based on Purser et al. (2003a) and describes details on how to ﬁnd
the RF coeﬃcients especially in the case of 4th-RF.

B.2 nth-order recursive ﬁlter
Deﬁne a diﬀerential operator D(n) as
σ 2 δx2 d2
1
D(n) = 1 −
+
2
2 dx
2!

Ã

σ 2 δx2 d2
2 dx2

!2

Ã

σ 2 δx2 d2
1−+ ··· +
n!
2 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(∞) givesÃσ 2 δx2 d2
D(∞) = exp −
2 dx2

!

(B.5)

−1
Let A(x) and C(x) be input and output signals of D(∞), respectively. Consider theirspectral representations Ã(k) and C̃(k):Z

A(x) =
C(x) =

1
Ã(k)eikx dk
2πZ
1C̃(k)eikx dk
2π

(B.6)
(B.7)

Then,
D(∞) C(x) =

1
2π

=

1
2π

Z
Z

Ã

σ 2 δx2 d2
C̃(k) exp −
2 dx2
Ã

!

eikx dk(B.8)!

σ 2 δx2 2 ikx
C̃(k) exp
k e dk2(B.9)

Thus, the spectral representation of the pseudo-diﬀerential operator D(∞) becomesÃσ 2 δx2 2kD̃(∞) = exp
2!

(B.10)For an input signal with the form of delta function
A(x) = δ(x) =

1
2π

184

Z

eikx dk(B.11)

−1
the transformation by D̃(∞)
is written as

−1
D̃(∞)
δ(x)

Z

Ã

!

σ 2 δx2 2 ikx
exp −k e dk
2

=

1
2π

=

1
x2
exp − 2 2
2π
2σ δx

Ã

(B.12)!

(B.13)

which is nothing but a Gaussian function with the length scale σδx. Thus, the trans−1−1formation D̃(∞)
gives a Gaussian ﬁlter, and D̃(n)is the nth-order approximation, thatis,−1
C = D(n)A(B.14)

gives an nth-order quasi-Gaussian ﬁlter that transforms the input signal A to get thesmoothed output C with the Gaussian shape with the length scale σδx.
Deﬁne a ﬁnite diﬀerence 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) inPurser et al. 2003a). Second order derivative with ﬁnite diﬀerencing is written asψi−1 − 2ψi + ψi+1
δx2
K
= − 2 ψi
δx

d2
ψi =
dx2

(B.16)(B.17)

where uniform grid spacing is assumed. Using eq.(B.17), we can rewrite eq.(B.4) as
σ2K
1
D(n) = 1 +
+
2
2!

Ã

σ2K
2

!2

1
+ ··· +
n!

Ã

σ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 K j (j = 1, · · · , n) give better results. In either185

way, D(n) is written in the form of an nth-degree polynomial w.r.t. K. The polynomial
can be written in a diﬀerent form asnYD(n) =

Ã

K
1−
κj

j=1

!

(B.19)

cf. eq.(A.1) in Purser et al. (2003a).
Deﬁne a shift operator Z asZψ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
K
1−
=
κj

Ã

1 − ζj Z −1
1 − ζj!Ã

1 − ζj Z
1 − ζj!

(B.22)

cf. eq.(A.5) in Purser et al. (2003a), where ζj satisﬁes
ζj2 − (2 − κj )ζj + 1 = 0(B.23)

Thus, the quasi-Gaussian transformation (eq.(B.14)) is rewritten as
Ci =

n
Y

Ã

j=1

1 − ζj1 − ζj Z −1!Ã

1 − ζj1 − ζj Z!

Ai(B.24)

We can separate the transformation process (eq.(B.24)) as follows:
Bi =
Ci =

n
Y

Ã

j=1
Ã
n
Y
j=1

1 − ζj1 − ζj Z −1
1 − ζj1 − ζj Z!

Ai(B.25)!

Bi(B.26)

αj Z −j Bi(B.27)

αj Z j Ci(B.28)

which are equivalent to
Bi = βAi +Ci = βBi +

n
X
j=1
n
X
j=1

186

whereα1 = ζ1 + ζ2 + · · · + ζn

(B.29)

α2 = ζ1 ζ2 + ζ1 ζ3 + · · · + ζn−1 ζn(B.30)

α3 = ζ1 ζ2 ζ3 + ζ1 ζ2 ζ4 + · · · + ζn−2 ζn−1 ζn(B.31)..
.
αn = ζ1 ζ2 ζ3 · · · ζn

(B.32)

β = (1 − ζ1 )(1 − ζ2 ) · · · (1 − ζn )= 1 − α1 − α2 − · · · − αn(B.33)

Eqs. (B.27) and (B.28) are nothing but forward and backward processes of the nth-RF(eqs.(refeq:nrﬀ2) 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 valuesof the coeﬃcients according to the length scale σ. The other polynomial (B.19) gives thecoeﬃcients
n
Y
j=1

Ã

K
1−
κj

!

µ

=1 +
µ

+

¶

1
1
1
− −
− ··· −
K
κ1 κ2
κn
¶

1
1
+ ··· +
K2
κ1 κ2
κn−1 κn..
.
µ+ (−1)n¶1Knκ1 κ2 · · · κn

(B.34)In the nth-RF, n roots of an nth-degree equation give κj . The details to get κj aredescribed in the next section in the case of the 4th-RF.

187

B.3 4th-order recursive ﬁlter
Since 4th-RF has been implemented in the 3DVAR system on the SPEEDY model, thissection describes speciﬁc 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 asD(4) = 1 +

σ2
K+
2

Ã

+

Ã

σ4 σ2
+
8
24

!

Ã

K2 +

!

σ6 σ4
σ2
K3
+
+
48 48 180
!

σ8
σ6
σ4
σ2
K4
+
+
+
384 192 1920 1120

(B.35)

Alternatively, eq.(B.34) in the 4th-order case gives
D(4) =

Ã

4
Y
j=1

K
1−
κj

!

µ

=1 +

¶

1
1
1
1
− −
−
−
K
κ1 κ2 κ3 κ4

µ

¶

1
1
1
1
1
1
+
+
+
+
+
+
K2
κ1 κ2 κ1 κ3 κ1 κ4 κ2 κ3 κ2 κ4 κ3 κ4
¶
µ
1
1
1
1
−
−
−
K3
+
−
κ1 κ2 κ3 κ1 κ2 κ4 κ1 κ3 κ4 κ2 κ3 κ4
µ
¶
1
+
K4
(B.36)
κ1 κ2 κ3 κ4
From the eqs.(B.35) and (B.36),
a =
b =
c =
d =

1
1
1
1
σ2
=− −
−
−
2
κ1 κ2 κ3 κ4
σ4 σ2
1
1
1
1
1
1
+
=
+
+
+
+
+
8
24
κ1 κ2 κ1 κ3 κ1 κ4 κ2 κ3 κ2 κ4 κ3 κ4
6
4
σ
σ
σ2
1
1
1
1
+
+
=−
−
−
−
48 48 180
κ1 κ2 κ3 κ1 κ2 κ4 κ1 κ3 κ4 κ2 κ3 κ4
σ8
σ6
σ4
σ2
1
+
+
+
=
384 192 1920 1120
κ1 κ2 κ3 κ4(B.37)(B.38)(B.39)(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ζj = 1 −
±
2

sµ

188

κj
1−
2

¶2

−1

(B.42)

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

α1 = ζ1 + ζ2 + ζ3 + ζ4

(B.43)

α2 = ζ1 ζ2 + ζ1 ζ3 + ζ1 ζ4 + ζ2 ζ3 + ζ2 ζ4 + ζ3 ζ4(B.44)

α3 = ζ1 ζ2 ζ3 + ζ1 ζ2 ζ4 + ζ1 ζ3 ζ4 + ζ2 ζ3 ζ4(B.45)

α4 = ζ1 ζ2 ζ3 ζ4(B.46)

β = 1 − α1 − α2 − α3 − α4(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: An Ensemble Adjustment Kalman Filter for Data Assimilation.
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. Meteorological Training Course Lecture Series, ECMWF, 75pp.
[8] Burgers, G., P. J. van Leeuwen, and G. Evensen, 1998: Analysis Scheme in theEnsemble 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 Paciﬁc sea level: Data assimilation via a reduced state space
Kalman ﬁlter. 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: Useof 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 thebred vectors young: Impact on data assimilation. The presentation ﬁle 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 Spectrumof 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 ﬁlter 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 andthree 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, Colorado, 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 ﬁltering. Tellus,
56A, 273-277.
[32] Ide, K., P. Courtier, M. Ghil, and A. C. Lorenc, 1997: Uniﬁed 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. Cambridge 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 Predictability, Reading, United Kingdom, ECMWF, 1-18.
[40] Lorenz, E. and K. Emanuel, 1998: Optimal sites for supplementary weather observations. J. Atmos. Sci., 55, 399-414.
[41] Marshall, J. and F. Molteni, 1993: Toward a Dynamical Understanding of PlanetaryScale Flow Regimes. J. Atmos. Sci., 50, 1792-1818.
[42] Miyoshi, T., 2004: Classical methods tour of advanced data assimilation using Lorenz
96 model. A ﬁnal report submitted to METO658E, Department of Meteorology, University of Maryland, College Park.
[43] Miyoshi, T., 2004: A review of the deterministic ensemble Kalman ﬁltering methods 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 andstochastically seeded bred vectors in the Lorenz 96 model with orography. In preparation.
[45] Molteni, F., 2003: Atmospheric simulations using a GCM with simpliﬁed physicalparameterizations. 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 ﬁltering. 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 ﬁlter 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 Aspectsof 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 ﬁlter: Perfect modelexperiments 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 Perturbed Observations. Mon. Wea. Rev., 130, 1913-1924.
[57] Wang, X. and C. H. Bishop, 2003: A Comparison of Breeding and Ensemble Transform 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 PositiveNegative 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