arXiv:1410.1233v6 [] 29 Oct 2017

EnKF-C user guide
version 1.63.0

Pavel Sakov
June 19, 2014 – October 31, 2017

Contents
Introduction

5

1 EnKF

6

1.1

Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2

EnKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3

EnKF analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.3.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.3.2

Some schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
ETKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
DEnKF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.3.3

Some numerical considerations . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.4

Localisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.5

Asynchronous DA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.6

EnOI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2 EnKF-C

17

2.1

Design considerations

2.2

The workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3

Starting up: example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4

Parameter files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4.1

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Main parameter file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Global analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1

2.4.2

Model parameter file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4.3

Grid parameter file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4.4

Observation types parameter file . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.4.5

Observation data parameter file . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.5

File name conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.6

prep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6.1

Observation types, products, instruments, batches, readers

. . . . . . . . . . 25

Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Products . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Instruments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Batches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Readers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.7

2.8

2.6.2

Asynchronous DA / FGAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.6.3

Superobing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

calc

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.7.1

Observation functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.7.2

Interpolation of ensemble transforms . . . . . . . . . . . . . . . . . . . . . . . 30

2.7.3

Adaptive moderation of observations . . . . . . . . . . . . . . . . . . . . . . . 30

2.7.4

Moderation of spread reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.7.5

Innovation statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.7.6

Impact of observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.7.7

Multiple model grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.7.8

“Multi-scale” localisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.8.1

2.9

Capping of inflation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

DA tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.10 Point logs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.11 Use of innovation statistics for model validation . . . . . . . . . . . . . . . . . . . . . 38
2

2.12 Bias correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.13 Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.14 System issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.14.1 Memory footprint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.14.2 Exit action . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.14.3 Dependencies and compilation issues . . . . . . . . . . . . . . . . . . . . . . . 40
2.15 Possible problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.15.1 calc becomes too slow after increasing localisation radius . . . . . . . . . . . 40
Acknowledgments

41

References

43

Abbreviations

44

Symbols

45

3

License
EnKF-C
Copyright (C) 2014 Pavel Sakov and Bureau of Meteorology
Redistribution and use of material from the package EnKF-C, with or without modification, are
permitted provided that the following conditions are met:
1. Redistributions of material must retain the above copyright notice, this list of conditions and
the following disclaimer. 2. The names of the authors may not be used to endorse or promote
products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE AUTHORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

4

Introduction
EnKF-C aims to provide a compact generic framework for off-line data assimilation (DA) into largescale layered geophysical models with the ensemble Kalman filter (EnKF). Here “light-weight” has
higher priority than “generic”; that is, the code is not designed to cover every virtual possibility for
the sake of it, but rather to be expandable in practical (from the author’s point of view) situations.
Following are its other main features:
- coded in C for GNU/Linux platform;
- model-agnostic;
- can conduct DA either in EnKF or ensemble optimal interpolation (EnOI) mode;
- permits multiple model grids;
- can handle rectangular or curvilinear horizontal grids, z or sigma vertical layers.
EnKF-C is available from https://github.com/sakov/enkf-c. This user guide is a part of the
EnKF-C package. It is also available from http://arxiv.org/abs/1410.1233.
The user guide has two main sections. Section 1 overviews the basics of the EnKF; section 2
provides technical description of EnKF-C.

5

Chapter 1

EnKF
1.1

Kalman filter

The Kalman filter (KF) is the underlying concept behind the EnKF. It is rather simple if formulated
as the recursive least squares.
Consider the global (in time) nonlinear minimisation problem
{xai }ki=1 = arg min Jk (x1 , . . . , xk ),

(1.1)

{xi }ki=1

Jk (x1 , . . . , xk ) = (x1 − xf1 )T (Pf1 )−1 (x1 − xf1 )
+

+

k
X

i=1
k
X
i=2

[yi − Hi (xi )]T (Ri )−1 [yi − Hi (xi )]
[xi − Mi (xi−1 )]T (Qi )−1 [xi − Mi (xi−1 )] .

(1.2)

Here {xai }ki=1 is a set of k state vectors that minimise the cost function (1.2); indices i = 1, . . . , k
correspond to a sequence of DA cycles, so that x1 is the estimated model state at the first cycle
and xk is the estimated model state at the last cycle; yi are observation vectors related to the
model state by nonlinear observation operators Hi (x); Mi (x) are nonlinear model functions that
map the model states xi−1 to xi ; Pf1 is the initial state error covariance; Ri are observation error
covariances; and Qi are model error covariances; and (·)T denotes matrix transposition.
The minimisation problem (1.1, 1.2) is, generally, very complicated, but, luckily, has an exact
solution in the linear case; moreover, this solution is recursive. Namely, assume that M and H are
affine:
Mi (x(1) ) − Mi (x(2) ) = Mi (x(1) − x(2) ),

(1.3a)

Hi (x(1) ) − Hi (x(2) ) = Hi (x(1) − x(2) ),

(1.3b)

where x(1) , x(2) are arbitrary model states, and Mi , Hi = Const. Then the cost function (1.2)
6

becomes quadratic and can be written in canonical form in regard to xk as follows:
Jk (x1 , . . . , xk ) = (xk − xa )T (Pa )−1 (xk − xa ) + J˜k−1 (x1 , . . . , xk−1 ),
k

k

k

so that
min Jk (x1 , . . . , xk ) = (xk − xak )T (Pak )−1 (xk − xak ) + Const.

{xi }k−1
i=1

(1.4)

(Proposition) Then
min Jk+1 (x1 , . . . , xk , xk+1 ) = (xk+1 − xak+1 )T (Pak+1 )−1 (xk+1 − xak+1 ) + Const,

{xi }ki=1

(1.5)

where
h
i
xak+1 = xfk+1 + Kk+1 yk+1 − Hk+1 (xfk+1 ) ,
Pak+1 = (I − Kk+1 Hk+1 )Pfk+1 ,

(1.6a)
(1.6b)

where

h
i−1
Kk+1 ≡ Pfk+1 (Hk+1 )T Hk+1 Pfk+1 (Hk+1 )T + Rk+1
and
xfk+1 = Mk+1 (xak ),

Pfk+1 = Mk+1 Pak (Mk+1 )T + Qk+1 .

(1.6c)
(1.7a)
(1.7b)

This solution is known as the Kalman filter (KF, Kalman, 1960). Equations (1.7) describe advancing
the system in time and represent the stage commonly called “forecast”, while equations (1.6)
describe assimilation of observations and represent the stage called “analysis”. The superscripts
f and a are used hereafter to refer to the forecast and analysis variables, correspondingly. The
forecast and analysis model state estimates xf and xa are commonly called (simply) forecast and
analysis. Matrix K is called Kalman gain.
There are a few things to be noted about the KF:
1. It follows from the Kalman filter (equations 1.4-1.7) that the state of the DA system (SDAS)
X is carried by the estimated model state vector and model state error covariance:
Xk = {xk , Pk }.

(1.8)

2. The KF provides solution for the last analysis, corresponding to xak in (1.1) (or, with a minor
re-formulation, to the last forecast); finding the full (global in time) solution requires application of the Kalman smoother (KS). Both the KF and KS can be derived by decomposition
of the positive (semi)definite quadratic function (1.2).
3. Because the SDAS represents a (part of a) solution of the global least squares problem, it
does not depend on the order in which observations are assimilated or on their grouping.
4. Ditto, the SDAS does not depend on a linear non-singular transform of the model state in the
sense that the forward and inverse transforms commute with the evolution of the DA system.
5. Solution (1.7, 1.6) can be used in a nonlinear case by approximating
Mi ← ∇Mi (xai−1 ),
Hi ← ∇Hi (xfi ),

in which case it is called the extended Kalman filter (EKF).
7

1.2

EnKF

The standard form of the KF (1.7, 1.6) is not necessarily the most convenient or suitable one in
practice. The corresponding algorithms can be prone to loosing the positive definiteness of the
state error covariance P due to round-up errors; and more importantly, explicit use of P makes
these algorithms non-scalable in regard to the model state dimension.
Both these immediate problems can be addressed with the ensemble Kalman filter, or the EnKF. In
the EnKF the SDAS is carried by an ensemble of model states E, which can be split into ensemble
mean and ensemble anomalies:
X = {E} = {x, A}.

(1.9)

It is related to the SDAS of the KF (1.8) as follows:
1
E 1,
m
1
P=
AAT ,
m−1
A ≡ E − x1T ,

(1.10a)

x=

(1.10b)
(1.10c)

where 1 is a vector with all elements equal to 1. The above means that the model state estimate
is given by the ensemble mean, while the model state error covariance P is implicitly represented
by the ensemble anomalies A via the factorisation (1.10b).
Representing the state error covariance via ensemble anomalies yields a number of numerical benefits. In large-scale geophysical systems the state size (∼ 105 − 109 ) makes it impossible to store and
manipulate the state error covariance P directly. At the same time it is often/typically possible
to represent essential variability via an ensemble of much smaller size (∼ 102 ) and manipulate P
implicitly via operations with A. Further, using A ensures positive semidefiniteness of P.
Storing the SDAS via an ensemble of model states is the first essential feature of the EnKF. In
theory, one could use it in an implementation of the KF along with explicitly calculated Jacobians
M and H in equations (1.7) and (1.6). The EnKF makes a further step and uses the ensemble
form of the SDAS for a derivative-less formulation of the KF. Moreover, it approximates derivatives
using ensemble of finite spread that characterises the estimated uncertainty in the state:
Hx
HA
Mx
MA
HMA

→

H(E) 1/m

(1.11a)
T



→

H(E) I − 11 /m

→

M(E) I − 11T /m

→

→

M(E) 1/m



T

(1.11b)
(1.11c)
(1.11d)


H ◦ M(E) I − 11 /m .

(1.11e)

This formulation is not the only one possible; one might also use the finite difference approximations:

8

Hx
HA
Mx
MA
HMA

→

→

→

→

→

H(x1T + εA) 1/m
T

(1.12a)
T



H(x1 + εA) I − 11 /m /ε

M(x1T + εA) 1/m
T

T



M(x1 + εA) I − 11 /m /ε

(1.12b)
(1.12c)
(1.12d)

H ◦ M(x1T + εA) I − 11T /m /ε.


(1.12e)

In this form the filter represents a derivative-less ensemble formulation of the EKF. In practice the
difference between the EnKF formulation (1.11) and the EKF formulation (1.12) is that the latter
is more sensitive to small-scale variability and therefore more prone to instability, similar to the
difference in behaviour of the Newton and secant methods.
The forecast stage of the EnKF involves only propagation of each ensemble member:
Efi = Mi (Eai−1 ).

(1.13)

This is a remarkably simple equation compared to the KF forecast equations (1.7), even though the
model error still needs to be handled in some way. Because propagation of each ensemble member
is independent from other members, the forecast stage in the EnKF is naturally parallelisable.
One way to handle model error in the EnKF is to include stochastic model error into the model
operator in (1.13). (This would make it different to the model operator in the KF, which is
deterministic.) Another option is to use the multiplicative inflation. The third option is to mimic
the treatment of model error in the KF, although this may require the “rank reduction” step
(Verlaan and Heemink, 1997, eq. 28).
At the analysis stage one has to update the ensemble mean and ensemble anomalies to match (1.6).
This involves handling ensemble as a whole, which is different to the forecast stage, when each
ensemble member is propagated individually.
Note that factorisation (1.10b) is not unique: if A satisfies (1.10b), then Ã = AU, where U is an
arbitrary orthonormal matrix UUT = I, also satisfies (1.10b). However, Ã should not only factorise
P, but also remain an ensemble anomalies matrix, Ã1 = 0. This requires an additional constraint
U1 = 1. Summarising, if E = x1T + A is an ensemble that satisfies (1.10), then ensemble
Ẽ = x1T + AUp ,

Up : Up (Up )T = I, Up 1 = 1

(1.14)

also satisfies (1.10). If E is full rank (i.e. rank(E) = min(m, n), where n is the state dimension), then
each unique Up generates a unique ensemble, and (1.14) describes all possible ensembles matching
a given SDAS of the KF. Such transformation of the ensemble is called ensemble redrawing. In the
linear case (i.e. for affine model and observation operators) redrawing of the ensemble in the EnKF
does not affect evolution of the underlying KF; and conversely, in the nonlinear case the redrawing
does indeed affect evolution of the underlying KF.

1.3

EnKF analysis

In this section we will give a brief overview of solutions for the EnKF analysis, and then describe
the particular schemes used in EnKF-C.
9

1.3.1

Overview

In the “baseline” EnKF (full-rank ensemble, no localisation) the analysed SDAS matches that of the
KF, although the algebraic side is indeed different. The update of the ensemble mean is generally
straightforward, in accordance with that in the KF (1.6a). The details may depend on the chosen
algorithm to achieve better numerical efficiency (see sec. 1.3.3).
The update of ensemble anomalies can be done either via a right-multiplied or left-multiplied (or
post/pre-multiplied) transform of the ensemble anomalies:
Aa = TL Af ,

(1.15)

Aa = Af TR .

(1.16)

or

TL and TR are referred to hereafter as left-multiplied and right-multiplied ensemble transform
matrices (ETMs), respectively. Note that for a full rank ensemble TR has to satisfy TR 1 = 1. It
follows from (1.14) that if TR is a particular solution for the right-multiplied ETM, then (for a full
rank ensemble) any other solution can be written as
T̃R = TR Up ,

Up : Up (Up )T = I, Up 1 = 1.

(1.17)

Similarly, the analysis increment can be represented as a linear combination of the forecast ensemble
anomalies:
xa = xf + Af w.

(1.18)

Equations (1.16) and (1.18) can be combined into a single transform of the ensemble:
Ea = Ef X5 ,
X5 =

(1.19)


1 T
1
11 + I − 11T
m
m



w1T + TR =






1
1 T
11 + w1T + I − 11T TR ,
m
m

(1.20)

as 1T w = 0. Because TR 1 = 1, for symmetric transforms (including ETKF and DEnKF) 1T TR =
1T , and (1.20) becomes
X5 = w1T + TR .

(1.21)

The designation X5 is used for historic reasons, following Evensen (2003).

1.3.2

Some schemes

As follows from the previous section, there are multiple solutions for the ETM that match the KF
covariance update equation (1.6b); however the particular solutions may have different properties
in practice due to the DAS nonlinearity, their algorithmic convenience, or their robustness in
suboptimal conditions. This section provides some background for the schemes used in EnKF-C:
- ETKF;
- DEnKF.
10

ETKF
It is easy to show using the definition of K (1.6c) and matrix shift lemma (1.23) that
(I − KH) Pf = (I − KH)1/2 Pf (I − KH)T/2 ,
which yields the following solution for the left-multiplied ETM:
TL = (I − KH)1/2

(1.22)

(Sakov and Oke, 2008b), that is
Aa = (I − KH)1/2 Af .

(1.22a)

Hereafter by X1/2 we denote the unique positive definite square root from a positive definite (generally, non-symmetric) matrix X, defined as X1/2 = VL1/2 V−1 , where X = VLV−1 is the eigenvalue
decomposition of X. By “matrix shift lemma” we refer to the following identity:
F(AB) A = A F(BA),

(1.23)

where F is an arbitrary function expandable into Taylor series. Rewriting (1.22a) as

1/2
1
a
f
f T
f T
−1
Af
A = I−
A (HA ) (HP H + R) H
m−1
and using the matrix shift lemma, we obtain:
1/2

1
f T
f T
−1
f
a
f
(HA ) (HP H + R) HA
A =A I−
m−1
which yields the corresponding to (1.22) right-multiplied ETM:

1/2
1
f T
f T
−1
f
TR = I −
(HA ) (HP H + R) HA
m−1

(1.24)

(Evensen, 2004). Applying the matrix inversion lemma
(A + ULV)−1 = A−1 − A−1 U(L−1 + VA−1 U)−1 VA−1 ,

(1.25)

(1.22) can be transformed to:
TL = (I + Pf HT R−1 H)−1/2

(1.26)

(Sakov and Bertino, 2011); and applying the matrix shift lemma yields the corresponding rightmultiplied ETM:
−1/2

1
f T −1
f
(HA ) R HA
,
(1.27)
TR = I +
m−1
also known as the ensemble transform Kalman filter, or ETKF (Bishop et al., 2001).
Historic reference. Another (and probably the first) solution for TR equivalent to (1.24) and (1.27) was found by Andrews
(1968):

−1
1
(HAf )T M−1/2 M1/2 + R1/2
HAf ,
(1.28)
TR = I −
m−1

11

where M ≡ HPf HT + R.

Equations (1.24, 1.27, 1.28) yield algebraically different expressions for the (unique) symmetric
right-multiplied solution. Apart from being the only symmetric solution, it also represents the
minimum distance solution for the ensemble anomalies: its ensemble of analysed anomalies is
closer to the ensemble of forecast anomalies with the inverse forecast (or analysis) covariance as the
metric than any other ensemble of analysed anomalies given by (1.17) (Ott et al., 2003, rev. 2005).
This means that in the above sense the symmetric right-multiplied solution preserves the identities
of ensemble members during analysis in the best possible way.
Note that while the left-multiplied solutions (1.22, 1.26) correspond to the symmetric right-multiplied
solution, they are not symmetric.
In a typical DAS with a large scale model one can expect m = 100, p = 103 − 107 , n = 106 − 109 ;
that is
m ≪ p ≪ n.

(1.29)

Therefore, considering the size of ETMs (n × n for left-multiplied ETMs and m × m for rightmultiplied ETMs), only right-multiplied solutions are suitable for use with large scale models. The
ETKF solution (1.27) represents the most popular option due to its simple form and numerical
effectiveness: for a diagonal R, it only requires to calculate inverse square root of a symmetric
m × m matrix. Also, along with the left-multiplied solution (1.26), it generally has better numerical
properties than solutions (1.22) and (1.24) due to the fact that the inverse square root in it is
calculated from the sum of a positive definite and a positive semi-definite matrices.

DEnKF
Assuming that KH is small in some sense, one can approximate solution (1.22) by expanding it
into Taylor series about I and keeping the first two terms of the expansion:
1
TL = I − KH.
2

(1.30)

This approximation is known as the deterministic Kalman filter, or DEnKF (Sakov and Oke,
2008a). It has a simple interpretation of using half of the Kalman gain for updating the ensemble anomalies; but apart from that the DEnKF often represents a good practical choice due
to its algorithmic convenience and good performance in suboptimal situations. The DEnKF is the
default scheme in EnKF-C.

1.3.3

Some numerical considerations

Instead of using the forecast ensemble observation anomalies HAf and innovation y − H(xf ) it is
convenient to use their standardised versions:
h
i √
s = R−1/2 y − H(xf ) / m − 1,
(1.31)
√
(1.32)
S = R−1/2 HAf / m − 1.
12

Then
w = Gs;

(1.33)

for the ETKF
TR = (I + ST S)−1/2 ,

(1.34)

and for the DEnKF
1
TR = I − GS,
2

(1.35)

where
G ≡ (I + ST S)−1 ST ,
T

T −1

= S (I + SS )

(1.36)

.

(1.37)

Here (1.36) involves inversion of an m × m matrix, while (1.37) involves inversion of a p × p matrix.
Therefore, in the DEnKF it is possible to calculate w and T using a single inversion of either a p × p
or m × m matrix, depending on the relation between the number of observations and the ensemble
size. In contrast, the ETKF (1.34) requires calculation of the inverse square root of an m × m
matrix. Then, one can use expression (1.36) for G and calculate both inversion in it and inverse
square root in (1.34) from the same singular value decomposition (SVD). This makes the DEnKF
somewhat more numerically effective because, firstly, one can exploit situations when p < m to
invert a matrix of lower dimension and, secondly, it requires only matrix inversion, which can be
done via Cholesky decomposition instead of SVD.

1.4

Localisation

Localisation is a necessary attribute of the EnKF systems with large-scale models, aimed at overcoming the rank deficiency of the ensemble. It can also be seen as aimed at reducing spurious long
range correlations occurring due to the finite size of the ensemble; or at limiting the impact of
distant observations because of the unreliability of the corresponding covariances.
There are two common localisation methods for the EnKF – covariance localisation (CL, Hamill and Whitaker,
2001; Houtekamer and Mitchell, 2001), also known as covariance filtering, and local analysis (LA,
Evensen, 2003; Ott et al., 2003, rev. 2005). Although CL may have advantages in certain situations
(non-local observations, “strong” assimilation), in practice the two methods produce similar results
(Sakov and Bertino, 2011). For algorithmic reasons EnKF-C uses LA.
Instead of calculating the global ensemble transform X5 , LA involves calculating local ensemble
i
transforms X5 for each element i of the state vector. This is done using local normalised ensemble
i
i
observation anomalies S and local normalised innovation s, obtained by tapering global S and s:
i

i

s ≡ s ◦ f,
i

(1.38a)
i

S ≡ S ◦ (f 1T ),

(1.38b)

13

i

where f is the vector of taper coefficients for element i, and A◦B denotes by-element, or Hadamard,
or Schur product of matrices A and B. We consider non-adaptive localisation only, when a taper
i
coefficient is a function of locations of the state element i (denoted as r) and observation o (denoted
i

{o}

i {o}

as r): fo = g(r, r), where g is the taper function. In layered geophysical models g is often assumed
to depend only on horizontal distance between these locations:
i

i

{o}

fo = g(|ρ − ρ|),

(1.39)
i

i

{o}

i

{o}

or on combination of horizontal and vertical distances, e.g.: fo = gxy (|ρ − ρ|)gz (|z − z|), where
r = (ρ, z), and ρ = (x, y). In the case (1.39) for a given set of observations the local ensemble
i
transform X5 depends only on horizontal grid coordinates of the state element xi and can be used
for updating all state elements with the same horizontal grid coordinates. This is currently the
only option in EnKF-C.
Smooth taper functions have advantage over non-smooth functions (such as the boxcar, or step
function) because they maintain the spatial continuity of the analysis. EnKF-C uses the popular
polynomial taper function by Gaspari and Cohn (1999), which has a number of nice properties.

1.5

Asynchronous DA

Observations assimilated at each cycle in the KF are assumed to be made simultaneously at the
time of assimilation. In such cases observations and DA method are referred to as synchronous. In
reality, observations assimilated at a given cycle are made over some period of time called “data
assimilation window” (DAW). If the DA method accounts for the time of observations, observations
and DA method are referred to as asynchronous.
The EnKF can be naturally extended for asynchronous DA. Let us consider the minimisation
problem (1.1, 1.2) in the case of perfect model Q = 0. It becomes
xa1 = arg min J(x1 ),
J(x1 ) = (x1 − xf1 )T (Pf1 )−1 (x1 − xf1 ) +
xi+1 = M(xi ),

i = 1, . . . , k − 1.

(1.40)
k
X
i=1

[yi − Hi (xi )]T (Ri )−1 [yi − Hi (xi )] ,

(1.41)
(1.42)

Compared to the original problem, the dimensionality of the solution is much reduced due to
relations (1.42), which mean that the model state at any time can be found by propagating the
initial state: x2 = M2 (x1 ), x3 = M3 ◦ M2 (x1 ), . . . . The cost function (1.41) can then be written
as
J(x1 ) = (x1 − xf1 )T (Pf )−1 (x1 − xf1 ) + [y − H ◦ M(x1 )]T R−1 [y − H ◦ M(x1 )] ,

(1.43)

where observations y represent the augmented observation vector: y = [y1T , . . . , ykT ]T , R is the
corresponding observation error covariance, and forward operator H ◦ M(x1 ) relates the initial
state x1 to observations.
Note that usually only H is a function that depends on assimilated observations; now, if we interpret
H ◦ M(x) = H[M(x)], we have to assume that M also depends on observations, propagating the
14

initial state to the time of each observation. It is also possible to interpret M(x) as the trajectory
starting from x, while H maps it to observations.
Apart from the operator H ◦ M, the cost function (1.43) has the same form as that for a single DA
cycle with synchronous observations. Consequently, in the linear case (1.3) one can use solutions
for w and T from section 1.3.3, subject to extending definitions of s (1.31) and S (1.32) as follows:
h
i √
(1.44)
s = R−1/2 y − H ◦ M(xf1 ) / m − 1,
√
(1.45)
S = R−1/2 H ◦ M Af1 / m − 1,
where H ◦ M is the tangent linear operator of H ◦ M about xf1 . This means that to account for
the time of observations in the EnKF one simply needs calculate innovation and forecast ensemble observation anomalies using ensemble at the time of each observation. There are no specific
restrictions on R, so that in theory observation errors can be correlated in time.
The minimisation problem (1.40),(1.43) implies assimilation time t = t1 ; however, in the linear
case the standardised innovation and ensemble observation anomalies in form (1.44,1.45) represent
objects invariant to assimilation time: the reference to x1 is only needed to define forward operator H ◦ M. Consequently, the ensemble transform X5 calculated from s and S can be applied to
ensemble at any particular time to yield (the same) analysed trajectories for the ensemble members: M(E) X5 = M(E X5 ). This time invariance of ensemble transforms can be used to update
ensemble back in time using observations from future cycles without the need in backward model
(Evensen and van Leeuwen 2000, sec. 6, Evensen 2003, app. D).
Note. The background term (x1 − xf1 )T (Pf )−1 (x1 − xf1 ) in (1.43) can be seen as accumulating the previous history of the
system rather than characterising the initial uncertainty in the global problem. In this case it is natural to anchor it to the
previous analysis:
J(x) = (x − xf )T (Pf )−1 (x − xf ) + [y − H ◦ M(x)]T R−1 [y − H ◦ M(x)] .
(1.46)
Here xf is the analysed state at the end of the previous cycle, considered as the forecast for the current cycle, and Pf is the
corresponding state error covariance. Minimising J yields the analysed initial model state xa , which in turn yields the analysed
trajectory; the analysed state error covariance is defined so that the analysed background term absorbs the observation term:
xa = arg min J(x),
Pa : J(xa ) = (xa − xf )T (Pa )−1 (xa − xf ) + Const.
This framework is a natural extension of the problem (1.1) in the linear, perfect-model case to continuous time; and is convenient
for iterative minimisation.

1.6

EnOI

The EnOI, or ensemble optimal interpolation (Evensen, 2003), can be defined as the EnKF with a
static or, more generally, pre-defined, ensemble anomalies. It can be summarised as follows:
xbi = Mi (xai−1 ),

xai

=

xbi

(1.47)

b

+ A wi ,

(1.48)

where xb is the forecast model state estimate referred to as background, and Ab is an ensemble of
static, or background, anomalies; the corresponding state error covariance Pb is also often referred
to as the background covariance.
The main incentive for using the EnOI is its low computational cost due to the integration of
only one instance of the model. Despite of the similarity with the EnKF, the EnOI is a rather
15

different concept, as there is no global in time cost function associated with it. Conceptually the
EnOI is closer to 3D-Var, as both methods use static (anisotropic, multivariate) covariance. It
is an improvement on the optimal interpolation, which typically uses isotropic, homogeneous and
univariate covariance.
In contrast to the EnKF, due to the use of a static ensemble the EnOI avoids potential problems
related to the ensemble spread; but at the same time it does critically depend on the ensemble,
while the EnKF with a stochastic model typically “forgets” the initial ensemble over time.
The EnOI can account for the time of observations by calculating innovation using forecast at
observation time, as in (1.44). This approach is commonly known as “first guess at appropriate
time”, or FGAT.

16

Chapter 2

EnKF-C
2.1

Design considerations

EnKF-C is designed to use horizontal localisation only. While some may argue that using vertical
localisation might help to decrease the ensemble size, we believe that, for example, in the ocean the
vertical structure is too complicated and non-uniform to allow simple and robust solutions in this
regard. Generally, dynamical processes include a variety of barotropic and baroclinic components,
and introducing vertical localisation in one form or another can be detrimental for the model’s
balances. On the other hand, with an ensemble size of about 100, normally one can ignore the
problem of spurious vertical correlations and leave the system to deal with the vertical covariances
on its own.
With the system using horizontal localisation only the model state effectively becomes a collection
of independent horizontal fields updated based on their correlations with local observations. The
assimilation is conducted by calculating a common horizontal field of local ensemble transforms
and applying them to each horizontal field of the model.

2.2

The workflow

EnKF-C conducts data assimilation in three stages: prep, calc and update.
prep preprocesses observations so that they are ready for DA. It has the following work flow:
- read original observations and fill each observation into a uniform structure measurement;
- write observations to file observations-orig.nc;
- combine observations into superobservations;
- write super observations to observations.nc.
calc calculates ensemble transforms for updating the forecast ensemble of model states (EnKF) or
the background model state (EnOI) in the following stages:

17

- read observations from observations.nc;
- calculate ensemble of forecast observations HEf (EnKF) or ensemble of background observation anomalies HAf and background observations Hxf (EnOI);
- for each horizontal grid cell calculate local ensemble transforms X5 (EnKF) or background
update coefficients w (EnOI);
- save these transforms to X5.nc (EnKF) or w.nc (EnOI);
- calculate and report forecast and analysis innovation statistics;
- calculate observation impact metrics DFS and SRF (sec. 2.7.6) and save them to enkf_diag.nc;
- save ensemble, observations and transforms/weights for specified horizontal locations to pointlog files (sec. 2.10).
Apart from this main mode of operation, calc can also calculate forecast innovations only or operate
in a single observation experiment mode.
update updates the ensemble (EnKF) or the background (EnOI) using transforms calculated by
calc, and saves the analysis according to specifications. It also adds analysis information to the
pointlog files, and can calculate and save the ensemble spread.

2.3

Starting up: example 1

It may be a good idea to start getting familiar with the system by running the example in
examples/1. The example has been put up based on the runs of regional EnKF and EnOI reanalysis systems for Tasman Sea developed by Bureau of Meteorology. It allows one to conduct
a single assimilation for 23 December 2007 (day 6565 since 1 January 1990) with either EnKF or
EnOI. To reduce the size of the system, the model state has been stripped down to two vertical
levels and 100 × 100 horizontal grid. Due to its size (almost 80 MB) the data for this example is
available for download separately from the EnKF-C code – see examples/1/README for details.

2.4

Parameter files

EnKF-C requires 5 parameter files to run:
- main parameter file;
- model parameter file;
- grid parameter file;
- observation types parameter file;
- and observation parameter file.
Examples of these parameter files can be found in examples/1. Running EnKF-C binaries with
--describe-prm-format in the command line provides information on the parameter file formats.

18

2.4.1

Main parameter file

The main parameter file specifies the main parameters of DA and 4 other parameter files. Its format
is described by running enkf_prep, enkf_calc or enkf_update with option --describe-prm-format:
>./bin/enkf_prep --describe-prm-format
Main parameter file format:
MODE
[ SCHEME
[ ALPHA
MODEL
GRID
OBSTYPES
OBS
DATE
[ WINDOWMIN
[ WINDOWMAX
ENSDIR
BGDIR
[ KFACTOR
[ RFACTOR
...
LOCRAD
[ STRIDE
[ SOBSTRIDE
[ FIELDBUFFERSIZE
[ INFLATION
...
[ ZSTATINTS
[ REGION
...
[ POINTLOG
...
[ EXITACTION
[ BADBATCHES
...

=
=
=
=
=
=
=
=
=
=
=
=
=
=

{ ENKF | ENOI }
{ DENKF* | ETKF } ]
<alpha> ]
(1*)
<model prm file>
<grid prm file>
<obs. types prm file>
<obs. data prm file>
<julian day of analysis>
<start of obs window in days from analysis> ] (-inf*)
<end of obs window in days from analysis> ]
(+inf*)
<ensemble directory>
<background directory>
(MODE = ENOI)
<kfactor> ]
(NaN*)
<rfactor> ]
(1*)

=
=
=
=
=

<loc. radius in km>
<stride> ]
<stride> ]
<fieldbuffersize> ]
<inflation> [ <VALUE>* | PLAIN ]

(1*)
(1*)
(1*)
(0.5*)

= [<z1> <z2>] ... ]
= <name> <lon1> <lon2> <lat1> <lat2> [[<z1> <z2>] ... ]
<i> <j> [grid name]]
= { BACKTRACE* | SEGFAULT } ]
= <obstype> <max. bias> <max. mad> <min # obs.> ]

Notes:
1. { ... | ... | ... } denotes the list of possible choices
2. [ ... ] denotes an optional input
3. ( ... ) is a note
4. * denotes the default value
5. < ... > denotes a description of an entry
6. ... denotes repeating the previous item an arbitrary number of times

Global analysis
It is possible to conduct global analysis by setting LOCRAD and STRIDE to large numbers. This is
demonstrated by target “global” in example 1.

19

2.4.2

Model parameter file

The model parameter file mainly describes the composition of the state vector by listing the model
variables and specifying the associated grids.
>./bin/enkf_prep --describe-prm-format model
Model parameter file format:
NAME

= <name>

VAR
[ GRID
[ INFLATION
[ RANDOMISE

= <name>
= <name> ]
(# grids > 1)
= <value> [<value> | PLAIN] ]
<deflation> <sigma> ]

[ <more of the above blocks> ]

Each model variable is described in a block started by the entry for the variable name. The inflation
parameters for a variable, if specified, override the common values set in the main parameter file
(sec. 2.8.1). Entry RANDOMISE makes it possible to specify a “forgetting” model for a variable
(sec. 2.12).
EnKF-C permits using multiple model grids. In this case each model variable must be associated
with one of the grids defined in the grid parameter file. See examples/4 for an example.

2.4.3

Grid parameter file

Grid parameter file describes grids used for model variables. Each grid is described in a section
started by the grid name entry and contains the grid name, grid data file, and names of the
dimensions and coordinates in the grid data file. It also contains variable names for the depth and
for number of layers in a vertical column (z grids) or land mask (sigma grids):
>./bin/enkf_prep --describe-prm-format grid
Grid parameter file format:
NAME
=
VTYPE
=
[ MAPTYPE
=
DATA
=
XDIMNAME
=
YDIMNAME
=
ZDIMNAME
=
XVARNAME
=
YVARNAME
=
ZVARNAME
=
DEPTHVARNAME
=
[ NUMLEVELSVARNAME
[ MASKVARNAME

<name> [ PREP | CALC ]
{ z | sigma }
{ binary* | kdtree } (curvilinear grids) ]
<data file name>
<x dimension name>
<y dimension name>
<z dimension name>
<x variable name>
<y variable name>
<z variable name>
<depth variable name>
= <# of levels variable name> (z) ]
= <land mask variable name> (sigma) ]

20

[ STRIDE
[ SFACTOR

= <stride> (common*) ]
= <spread factor> (1.0*) ]

[ <more of the above blocks> ]

The code is supposed to automatically identify the type of horizontal grid used, while the type of
the vertical grid has to be specified explicitly.
At the moment, EnKF-C supports 3 main types of horizontal grids:
- equidistant rectangular grids aligned with physical coordinates;
- non-equidistant rectangular grids aligned with physical coordinates;
- quadrilateral simply connected grids (via libgu).
For rectangular grids the code tries to determine and handle periodicity in X or Y directions.
With some effort EnKF-C can also handle topologically screwed grids (like tri-polar ORCA grids).
This is achieved by creating an auxiliary grid that is simply connected in grid index space by
masking some nodes (replacing their coordinates with NaNs) in the original grid. Because only
prep requires mapping from physical to grid space, this auxiliary grid is then used in prep only,
while the original grid is used in calc and update. To specify that a grid entry is to be used in prep
only one needs to add qualifier “PREP” after its name; and to specify that it is to be used in calc
and update one needs to add qualifier “CALC”. See examples/3 for an example.
From version 1.37.0 such grids can also be handled by using grid mapping via kd-tree instead of
the default mapping via spatial binary tree, which is now available from package gridutils-c. This
is done by setting parameter MAPTYPE to “k”.

2.4.4

Observation types parameter file

Observation types are the interface that connects model and observations. They are specified in
a separate parameter file. Each observation type is described in a separate section identified by
the entry NAME. Apart from the type name, the section must contain the tag for the observation
kind; associated model variable; the tag for the associated observation operator; and the allowed
range. The optional parameters include the R-factor for the type (sec. 2.9) and spatial limits on
corresponding observations. For example:
>./bin/enkf_prep --describe-prm-format obstypes
Observation types parameter file format:

[
[
[
[

NAME
ISSURFACE
VAR
OFFSET
MLD_VARNAME
MLD_THRESH
HFUNCTION
ASYNC

=
=
=
=
=
=
=
=

<name>
{0 | 1}
<model variable name> [...]
<file name> <variable name> ]
<model varname> ]
<threshold> ]
<H function tag>
<time interval> ]

(none*)
(none*)
(NaN*)
(synchronous*)

21

[
[
[
[
[
[
[
[
[
[
[

LOCRAD
WEIGHT
RFACTOR
MINVALUE
MAXVALUE
XMIN
XMAX
YMIN
YMAX
ZMIN
ZMAX

=
=
=
=
=
=
=
=
=
=
=

<locrad> ... ]
<weight> ... ]
<rfactor> ]
<minimal allowed
<maximal allowed
<minimal allowed
<maximal allowed
<minimal allowed
<maximal allowed
<minimal allowed
<maximal allowed

value> ]
value> ]
X coordinate>
X coordinate>
Y coordinate>
Y coordinate>
Z coordinate>
Z coordinate>

]
]
]
]
]
]

(global*)
(1*)
(1*)
(-inf*)
(+inf*)
(-inf*)
(+inf*)
(-inf*)
(+inf*)
(-inf*)
(+inf*)

[ <more of the above blocks> ]

The tags for available observation operators are listed in array allhentries in file calc/allhs.c.
The OFFSET entry may be used for adding the known model bias to observations, for example, to
specify the mean dynamic topography (MDT) when assimilating sea level anomaly (SLA) observations.
The localisation radius for an observation type, if specified, overrides the common value from
the main parameter file. The R-factors for each observation type are obtained by multiplying the
common value by the observation type value. (More on localisation radius and R-factor in sec. 2.9.)
The entry OFFSET allows one to add a new observation type without re-compiling the code. This
can be useful e.g. for evaluating a new instrument (by getting an independent statistics for the
new instrument and seeing the change in statistics for the old instruments).
The entries MLD_VARNAME and MLD_THRESH are used to calculate the model mixed layer depth for
projecting the surface bias.

2.4.5

Observation data parameter file

Observation data parameter file specifies observations to be assimilated. EnKF-C has a simple
policy in this regard: if a data file is listed in the observation data parameter file, then observations
from this file are assimilated. This allows one using custom observation windows for particular
observation types, instruments etc., specifying details on the script level during the parameter file
generation.
In practice some of observations specified in the observation data parameter file can be outside the
observation time window for the cycle. In this case the exact boundaries of the observation window
can be specified by entries WINDOWMIN and WINDOWMAX in the main parameter file; observations with
time outside interval [DATE-WINDOWMIN, DATE+WINDOWMAX) will be discarded.
The observation parameter file contains an arbitrary number of sections identified by entries
PRODUCT. Each section specifies the observation type, input files, reader and, possibly, observation error:
./bin/enkf_prep --describe-prm-format obsdata

22

Observation data parameter file format:
PRODUCT
READER
[ PARAMETER
...
TYPE
FILE
...
[ ERROR_STD
...

= <product>
= <reader>
<name> = <value> ]
= <observation type>
= <data file wildcard>
= { <value> | <data file> } [ EQ* | PL | MU | MI | MA ] ]

[ <more of the above blocks> ]

Observation files can be defined using wildcards “*” and “?”. Missing a file is reported in the log
and is not considered to be a fatal error.
The line in the above example starting with ERROR_STD specifies the observation error. It can
contain either a number or a file name. In the case of entering the file name there also should be
another entry in the same line specifying the name of the variable to be read. The variable should
have the same dimension (2D or 3D) as the associated observation kind as described by the field
issurface in the array allkinds in file common/obstypes.c (sec. 2.4.4).
The line with observation error can also have another token
p specifying the type of operation to be
2 + σ 2 ), MULT (σ
conducted: EQUAL (σtot ← σnow , default), PLUS (σtot ← σtot
tot ← σtot σnow ), MIN
now
(σtot = max(σtot , σnow )), or MAX (σtot = min(σtot , σnow )). There can be several error entries in a
section in the observation parameter file.
The observation time only matters if the observation type is specified to be “asynchronous” (see
sec. 2.6.2). In this case the model estimation for the observation is made by using model state at
the appropriate time. Otherwise, observations are assumed to be made at the time of assimilation,
regardless of the actual observation time.
Note that there can be multiple blocks with the same product. This enables custom treatment of
some specific data. For example, the following entries override observation error for Geosat (files
with prefix g1_) on 23 May 2006:
# set observation error for Geosat to 7cm
product == RADS
type = SLA
reader = standard2
file=/short/p93/pxs599/obs/RADS-IB/y2006/m05/g?_d23.nc
error_std = 0.07
# use default errors for other altimeters
product == RADS
type = SLA
reader = standard2
file=/short/p93/pxs599/obs/RADS-IB/y2006/m05/[!g]?_d23.nc

23

2.5

File name conventions

EnKF-C assumes that the ensemble and background file names have some predefined formats. The
file name for member mid and model variable varname is assumed to be sprintf("mem%03d_%s.nc",
mid, varname). The background file (EnOI only) for variable varname is assumed to be sprintf("bg_%s.nc",
varname). The above names are used for reading forecast states for synchronous DA and for
writing analyses, in the case if the analyses are appended to forecasts (true by default). For
asynchronous DA the member and background file names for the time slot t are assumed to
be sprintf("mem%03d_%s_%d.nc", mid, varname, t) and sprintf("bg_%s_%d.nc", varname,
t), correspondingly.

2.6

prep

prep is the first stage of data assimilation in EnKF-C. Its preprocesses observations by bringing
them to a common form and merging close observations into so called superobservations.
By design, prep is supposed to be light-weight, so that it does not read either the ensemble or
background, and the only model information it needs is the model grid. (Note that this may
require some additional processing at later stages for models with dynamic grid, such as HYCOM.)
The name of the binary (executable) for prep is enkf_prep. It has the following usage and options:
>./bin/enkf_prep
Usage: enkf_prep <prm file> [<options>]
Options:
--describe-prm-format [main|model|grid|obstypes|obsdata]
describe format of the parameter file and exit
--describe-superob <sob #>
print composition of this superobservation and exit
--consider-subgrid-variability
increase error for a superobservation when there is large subgrid variability
--log-all-obs
put all obs into observations.nc (default: obs within model domain only)
--no-superobing
--version
print version and exit

enkf_prep writes the preprocessed observations to file observatons.nc. It also writes original
(not superobed) observations to observatons-orig.nc (provided that superobing has not been
switched off by flag --no-superobing). By default, only observations within the model grid are
written to it, but the the option --log-all-obs changes this behaviour to writing all observations
from the input files.

24

2.6.1

Observation types, products, instruments, batches, readers

Types
Each observation has a number of attributes defined by the fields of the structure observation. One
of them is observation type, which characterises the observation in a general way and relates it to
the model state. For example, typical oceanographic observations may have tags SLA (for sea level
anomalies), SST (sea surface temperature), TEM (subsurface temperature) and SAL (subsurface
salinity). Different types can be related to the same model variable, as do SST and TEM in the
above example. Observation types are described in the corresponding parameter file (sec. 2.4.4).

Products
An observation is also characterised by “product”. It can be a tag for an organisation that provides
data from certain observational platforms, e.g.:
PRODUCT == RADS
TYPE = SLA
READER = standard2
FILE = obs/RADS-IB/y2007/m12/??_d19.nc
FILE = obs/RADS-IB/y2007/m12/??_d20.nc
FILE = obs/RADS-IB/y2007/m12/??_d21.nc
FILE = obs/RADS-IB/y2007/m12/??_d22.nc
FILE = obs/RADS-IB/y2007/m12/??_d23.nc
PRODUCT == NAVO
TYPE = SST
READER = standard
FILE = obs/NAVO/navo_20071219.nc
FILE = obs/NAVO/navo_20071220.nc
FILE = obs/NAVO/navo_20071221.nc
FILE = obs/NAVO/navo_20071222.nc
FILE = obs/NAVO/navo_20071223.nc

Instruments
The observational data from a product can be collected by a number of instruments. The corresponding field in the measurement structure is supposed to be filled by the observation reader.

Batches
An observation can be attributed to one of the groups called “batches”, such as altimeter passes,
Argo profiles etc., to enable detection and discarding of bad batches. Programmatically, to switch
on capabilities associated with observation batches for a particular kind of observations, the observation batch ID needs to be set by the corresponding observation reader.

25

A batch of observations is considered bad if either the magnitude of its innovation bias or the mean
magnitude of innovation exceed specified thresholds. Specifications for bad batches can be set in
the parameter file as follows:
BADBATCHES = SLA 0.06 0.10 500
BBADBATCHES = TEM 4 5 0
BADBATCHES = SST 0.5 2 10000
BADBATCHES = SAL 1.5 2 0

The above entry means that any batch of observations of type SLA (typically, an orbit) containing more than 500 observations and having either mean innovation greater than 0.06 (meter) in
magnitude or mean absolute innovation greater than 0.10 is considered to be bad. Similarly, a
TEM batch (typically, a profile) is considered bad if the mean innovation exceeds 4 (degrees) or the
mean absolute innovation exceeds 5 (degrees). The parameter file can have arbitrary number of
such entries. Information about bad batches is written by enkf_calc to the file badbatches.out.
When enkf_prep detects the presence of such file, it marks the corresponding observations as bad.
Therefore, the workflow for detecting and eliminating bad batches of observations is as follows:
1. specify bad batches in the parameter file;
2. make a pilot run of enkf_prep;
3. run enkf_calc with the flag --forecast-stats-only;
4. remove observations.nc and observations-orig.nc;
5. calculate analysis in a “normal” way by running enkf_prep, enkf_calc and enkf_update.
Note that bad batches are identified at stage 2 based on superobservations formed by observations
from the same batch only. Consequently, detection of bad batches may become unreliable if there
are not enough such “clean” superobservations. It is possible to alleviate this problem by (1)
switching off superobing by setting SOBSTRIDE = 0 or specifying --no-superobing for enkf_prep
and/or (2) by processing observations by parts, one type (or instrument) at a time, and repeating
stages 1 and 2 as many times as necessary.

Readers
The function of data readers is to read observations in specified files and parse them sequentially
into struct observation defined in common/observations.h.
At the moment there are few agreed standards for observational data. Even data of similar nature
may have specific features that require creating a new or modify an existing observation data reader.
It is therefore likely that developing a completely new system will require add new data readers.
The new readers can be created based on existing ones, or developed from scratch. The available
readers are listed by the variable allreaders defined in prep/allreaders.c. Each reader can
be specified in the observation data parameter file with an arbitrary number of parameters. For
26

example, the following section changes the default minimal depth for using altimetry observations
to 150 m:
(...)
product == RADS
type = SLA
reader = standard2
parameter MINDEPTH = 150
(...)
There are also several generic readers:
reader_xy_scattered
reader_xy_gridded
reader_xyz_scattered
reader_xyz_gridded
These readers allow one to specify a number of parameters by means of the entry PARAMETER of the
observation data parameter file. At the moment the only obligatory parameter for these readers is
VARNAME; a number of other obligatory parameters (e.g. LONNAME, LATNAME) can be specified either
via the entry PARAMETER or by using default variable names in data files.

2.6.2

Asynchronous DA / FGAT

An observation type can be specified as “asynchronous” by specifying entry ASYNC entry in the
observation types parameter file (sec. 2.4.4), e.g.:
NAME = SLA
(...)
ASYNC = 1
(...)

The above means that SLA and SST observations are considered to be asynchronous, with the
time quantification of 1 day. If, for example the assimilation time is specified as “6085.5 days since
1990-01-01”, then the SLA and SST observations will be binned into 1-day time intervals centred
at the time of assimilation, i.e. from day 6080.0 to day 6081.0, 6081.0 to 6082.0, and so on, and
estimated versus the corresponding model states.
The model states used to calculate forecast observations are matched by file names, which are
supposed to be of the form mem<xxx>_<variable name>_<time shift>.nc (for the EnKF) or
bg_<variable name>_<time shift>.nc (for the EnOI). Here “time shift” is the number of the
bin, with “0” corresponding to the bin centred at the time of assimilation, “-1” to the previous bin,
“1” to the next bin, and so on. If the corresponding members (or the background files, in the case
of EnOI) are found, the observations are assimilated asynchronously; if they are not found, then
the observations are assimilated synchronously. This can be tracked from the calc log file, e.g.:
27

calculating ensemble observations:
2014-03-22 06:28:28
ensemble size = 96
distributing iterations:
all processes get 6 iterations
process 0: 0 - 5
SST |aaaaaa|aaaaaa|aaaaaa|aaaaaa|aaaaaa
SLA |aaaaaa|aaaaaa|aaaaaa|aaaaaa|aaaaaa
TEM ......
SAL ......

The entries “a” mean that the observations are assimilated asynchronously. They would be replaced
by “s” if assimilated synchronously. The vertical lines indicate the time slots for asynchronous DA;
in the above example the DAW has 5 time slots. The entries “.” indicate calculating ensemble
observations for synchronous observations. Note that only the master process is writing to the
log here, which explains why there is only output from 6 members in the log above. For EnOI
there also will be “+” at the end of the line for each observation type indicating reading of the
corresponding background field.

2.6.3

Superobing

“Superobing” is the process of reduction of the number of observations by merging spatially close
observations before their assimilation. EnKF-C merges observations if:
- they belong to the same model grid cell;
- are of the same type;
- for asynchronous observations – belong to the same time slot.
The horizontal size of superobing cells can be increased from the default of 1 model grid cell to
N × N cells by setting SOBSTRIDE = <N> in the parameter file; the vertical size is always equal to
1 layer. Setting SOBSTRIDE = 0 switches superobing off.
The observations are merged by averaging their values, coordinates and times with weights inversely
proportional to the observation error variance. The observation error variance of a superobservation
is set to the inverse of the sum of inverse observation error variances of the merged observations.
The product and instrument fields of the superobservation are set either to those of the merged
observations or to -1, depending on whether the merged observations have the same values for these
fields or not.
Command line parameter --consider-subgrid-variability switches on considering the subgrid
variability by calculating standard deviation of the merged observations σsub and using σobs =
max(σobs , σsub ). The calculation of σsub is currently done in a rather crude way, assuming equal
weights for all merged observations.

28

2.7

calc

calc is the second stage of data assimilation in EnKF-C. It calculates 2D arrays of local ensemble
transforms X5 (for EnKF) or coefficients w (for EnOI).
The name of the binary for calc is enkf_calc. It has the following usage and options:
>./bin/enkf_calc
Usage: enkf_calc <prm file> [<options>]
Options:
--describe-prm-format [main|model|grid|obstypes]
describe format of a parameter file and exit
--forecast-stats-only
calculate and print forecast observation stats only
--ignore-no-obs
proceed even if there are no observations
--no-mean-update
update ensemble anomalies only
--point-logs-only
skip calculating transforms for the whole grid and observation stats
--print-batch-stats
calculate and print global biases for each batch of observations
--single-observation-xyz <lon> <lat> <depth> <type> <inn> <std>
assimilate single observation with these parameters
--single-observation-ijk <fi> <fj> <fk> <type> <inn> <std>
assimilate single observation with these parameters
--use-rmsd-for-obsstats
use RMSD instead of MAD when printing observation stats
--use-these-obs <obs file>
assimilate observations from this file; the file format must be compatible
with that of observations.nc produced by ‘enkf_prep’
--version
print version and exit

The option --forecast-stats-only can be used for quick calculation of the innovation statistics
for a given background (or ensemble). This can be used, for example, for obtaining the persistence
statistics, that is, the innovation statistics for the previous analysis.
The options --single-observation-xyz and --single-observation-ijk provide an easy way
to conduct the so called single observation experiments, with the observation coordinates provided
either in spatial or grid coordinates, correspondingly. Parameter <value> defines innovation rather
than the observation value. Normally, this experiments would be conducted in the EnOI mode,
calculating increment (option --output-increment of enkf_update) rather than analysis. When
run in the EnKF mode, the increment (or analysis, depending on specifications) for each member
is calculated.
Note that the calculated transforms do not incorporate inflation. Inflation is applied during update
according to specifications (sec. 2.8.1).

29

2.7.1

Observation functions

Model estimations for observations of each type are calculated using observation functions specified
for this type by entry HFUNCTIONS in the observation types parameter file, e.g.:
NAME = SLA
...
HFUNCTION = standard
...

The available functions for each observation type are specified by the variable allhentries in
calc/allhs.c. The “standard” functions do normally perform 2D or 3D linear interpolation from
the corner model grid nodes for the cell containing the observation.

2.7.2

Interpolation of ensemble transforms

Local ensemble transforms X5 (EnKF) or local ensemble weights w (EnOI) represent smooth fields
with the characteristic spatial variability scale of the localisation radius. This smoothness allows
one to reduce the computational load in calc by calculating local transforms or weights on a subgrid
with a specified stride only, and using linearly interpolated transforms or weights in the intermediate
grid cells. The value of the stride is defined by the STRIDE entry in the parameter file and can be
overwritten for a particular grid in the grid parameter file.

2.7.3

Adaptive moderation of observations

One of standard quality control (QC) procedures in DA is the so called background check, when an
observation is compared with the forecast and discarded if the innovation magnitude exceeds some
specified threshold. The downside of this approach is that it can not distinguish between situations
of an outlier, big model error (e.g. because of an error in forcing), or model divergence. While one
probably would like to discard an outlier, it is usually desirable to make use of valid observations,
although, perhaps, with a reduced impact, to avoid “over-stressing” the model. In EnKF-C this
is achieved by adaptive moderation of the observation impact by restricting the magnitude of the
increment from a given observation in observation space by K times magnitude of the spread of
the forecast ensemble.
Specifically, the adaptive moderation of the observation impact is conducted by smoothly increasing
the observation error depending on the magnitude of innovation as follows:
1/2

2 2
2
− σf2 ,
) + σf2 d2 /K 2
← (σf2 + σobs
σobs

where σobs is the observation error standard deviation, σf – forecast ensemble spread, d – innovation,
and K – the so called K-factor defined in the main parameter file (sec. 2.4.1). Tests with small
models show that setting the K-factor to say 2-3 or above makes a marginal impact (if any) on
performance of weakly suboptimal systems, while still can be quite beneficial in stressful situations.

30

2.7.4

Moderation of spread reduction

The moderating parameter α ∈ [0, 1] specified in the main parameter file via the entry ALPHA
allows one to reduce the contraction of ensemble during assimilation, while leaving the increment
unchanged. It modifies the DEnKF expression (1.35) for TR as
TR = I −

α
GS,
2

and the ETKF expression (1.34) as
TR = (I + α ST S)−1/2 .

2.7.5

Innovation statistics

In its course calc calculates some basic innovation statistics: number of observations, mean absolute
forecast innovation, mean absolute analysis innovation, mean forecast innovation, mean analysis
innovation, mean forecast ensemble spread, and mean analysis ensemble spread. This statistics is
provided for each region defined in the main parameter file (sec. 2.4.1), as well as for each time slot
defined for asynchronous DA, and for each instrument. By default, EnKF-C defines one statistical
region “Global” with domain [x1 , x2 ] = [−999, 999], [y1 , y2 ] = [−999, 999].
In addition, for 3D observations calc also calculates observation statistics in specified depth intervals. The global (common for all regions) intervals can be set by the entry ZSTATINTS in the main parameter file; by default, three intervals are defined: [0 DEPTH_SHALLOW], [DEPTH_SHALLOW DEPTH_DEEP],
and [DEPTH_DEEP DEPTH_MAX], where DEPTH_SHALLOW, DEPTH_DEEP and DEPTH_MAX are the macros
defined in common/definitions.h. These common intervals can be overriden by settings for particular statistical regions.
Following is an example of innovation statistics written to the log (standard output) of enkf_calc:
printing observation statistics:
region obs.type
# obs. |for.inn.| |an.inn.|
for.inn.
an.inn. for.spread an.spread
-----------------------------------------------------------------------------------------Tasman
SLA
3003
0.067
0.038
0.033
0.012
0.035
0.025
-4
712
0.058
0.038
0.035
0.013
0.028
0.021
-3
785
0.093
0.040
0.060
0.019
0.052
0.034
-2
700
0.062
0.043
0.030
0.016
0.027
0.021
-1
668
0.049
0.031
0.017
0.004
0.028
0.021
0
138
0.078
0.033
-0.043
-0.016
0.045
0.029
j1
1323
0.070
0.033
0.041
0.016
0.037
0.024
n1
876
0.073
0.042
0.052
0.025
0.036
0.026
g1
785
0.054
0.042
-0.004
-0.009
0.029
0.024
N/A
19
0.101
0.037
0.097
0.031
0.059
0.036
SST
9316
0.346
0.174
-0.215
-0.094
0.358
0.254
-4
2946
0.327
0.166
-0.236
-0.092
0.342
0.245
-3
2733
0.368
0.183
-0.270
-0.133
0.362
0.256
-2
2560
0.352
0.169
-0.167
-0.057
0.370
0.262
-1
580
0.342
0.191
-0.148
-0.093
0.414
0.291
0
497
0.305
0.182
-0.126
-0.075
0.307
0.225
AVHRR
9316
0.346
0.174
-0.215
-0.094
0.358
0.254

31

TEM
ARGO
0-50m
50-500m
>500m
SAL
ARGO
0-50m
50-500m
>500m

768
768
125
451
192
768
768
125
451
192

0.581
0.581
0.418
0.678
0.458
0.079
0.079
0.079
0.092
0.048

0.365
0.365
0.230
0.403
0.365
0.060
0.060
0.063
0.067
0.041

-0.245
-0.245
0.049
-0.266
-0.387
0.014
0.014
0.031
0.026
-0.027

-0.151
-0.151
0.027
-0.141
-0.291
0.019
0.019
0.035
0.032
-0.021

0.320
0.320
0.365
0.360
0.196
0.033
0.033
0.034
0.039
0.018

0.251
0.251
0.281
0.278
0.170
0.028
0.028
0.030
0.032
0.016

This excerpt shows innovation statistics for the region “Tasman”. It contains sections for SST,
SLA and TEM observations. The summary statistics for each observation type is shown at the
top of each section; then statistics for days -4, -3, -2, -1 and 0 of a 5-day DAW are shown for the
two asynchronous types, SST and SLA. After that, statistics for particular instruments is shown;
“N/A” corresponds to superobservations resulted from merging observations from two or more
instruments. For subsurface temperature also statistics for shallow (0–50 m), deep (>500 m), and
intermediate (50–500 m) observations is given.
The analysis innovation statistics is calculated from the updated (analysis) ensemble observations
by calc, thus avoiding the need to access analysis files produced later by update. The update of
ensemble observations is performed in the same way as that of any other element of the state vector:
for the EnKF – by applying the appropriate local ensemble transforms to the forecast ensemble
observations,
H(Ea ) ← H(Ef ) X5 ;
and for the EnOI – by applying the appropriate local linear combination of the ensemble observation
anomalies:
h
i
H(Ea ) ← H(xf ) + (HAf )w 1T + HAf .

2.7.6

Impact of observations

In the course of its work calc routinely calculates two metrics for assessing the impact of observations, degrees of freedom of signal (DFS) and spread reduction factor (SRF):
DFS = tr(KH) = tr(GS),
s
s
f
T
−1
tr(HP H R )
tr(ST S)
−1=
− 1,
SRF =
a
T
−1
tr(HP H R )
tr(GS)
where tr(·) is the trace function. The values of these metrics for each local analysis, calculated
both for all observations and for observations of each type only, are written to file enkf_diagn.nc.
Note that the in EnKF-C DFS and SRF are calculated from the above expressions and represent
theoretical values for the EnKF analysis; they coincide with the actual DFS and SRF values only
for the ETKF, but not for the DEnKF, which is an approximation of the KF (and indeed not for
the EnOI, which is not even an approximation).

32

In the EnKF context DFS is a useful indicator of potential rank problems. Normally, it should
not exceed a fraction (a half, or better, a quarter) of the ensemble size. SRF shows the “strength”
of DA. “Strong” DA implies a close to optimal system, which indeed never happens in practice.
Therefore, ideally, SRF should be small (below 1, on average).

2.7.7

Multiple model grids

EnKF-C permits using multiple model grids, in which case the ensemble transforms are calculated
sequentially for each of the grids. These transforms are then used for updating of the model
variables defined on the corresponding grids.

2.7.8

“Multi-scale” localisation

It is possible to specify the localisation taper function as a linear combination of the Gaspari and
Cohn’s taper functions with different support radii:
f (r) =

N
X

wi f0 (

i=1

r
),
Ri

where wi is the weight, r is the distance, Ri is the support radius, and

0 ≤ x ≤ 1,
 1 − 35 x2 + 58 x3 + 21 x4 − 14 x5 ,
x
2 −1
5 2
5 3
1 4
1 5
− x + 4 − 5x + 3 x + 8 x − 2 x + 12 x ,
1 < x ≤ 2,
f0 ( ) =
 3
2
0,
2 < x.

This can be set by entries LOCRAD and WEIGHT either in the main parameter file or in the observation
types parameter file, e.g.:
LOCRAD 150 500
WEIGHT 0.9 0.1
(recall that entries in the observation types parameter file for particular observation types override
the common settings in the main parameter file). Note that the weights are normalised so that
their sum is equal to 1.

2.8

update

update is the third and final stage of data assimilation in EnKF-C. It updates the ensemble (EnKF)
or the background (EnOI) by applying the transforms calculated by calc.
The name of the binary for update is enkf_update. It has the following usage and options:
>./bin/enkf_update
Usage: enkf_update <prm file> [<options>]

33

Options:
--calculate-spread
calculate ensemble spread and write to spread.nc
--describe-prm-format [main|model|grid]
describe format of a parameter file and exit
--direct-write
write fields directly to the output file (default: write to tiles first)
--joint-output
append analyses to forecast files (default: write to separate files)
--leave-tiles
do not delete tiles
--no-fields-write
do not write analysis fields
--output-increment
output analysis increment (default: output analysis)
--write-inflation
write adaptive inflation magnitudes to inflation.nc
--version
print version and exit

The option --joint-output tells update to append analyses to the corresponding forecast files,
using new variable names constructed by concatenating the forecast variable names and suffix _an.
By default the analyses are written to separate files, using the same variable names as the forecast
files, but with an extra suffix .analysis or .increment added to the file name, depending on
whether the analysis or increment is written.
By default, update first writes each updated horizontal field of the model to a separate file (referred
to here as a tile), and then concatenates these fields into analysis files. The tiles are removed
after writing the analysis files; one may save time for allocating them on disk in the next cycle by
leaving them on disk by using option --leave-tiles. This approach is somewhat less effective than
direct writing to analysis files (without intermediate tiles), but, unfortunately, the direct writing is
generally not reliable due to parallel I/O issues with NetCDF. Note that in some cases it proved
to be possible to obtain robust performance with direct write using “classic” or “64-bit-offset”
NetCDF formats.

2.8.1

Capping of inflation

Applying spatially uniform ensemble inflation involves areas with no local observations, where no
assimilation is conducted. It can gradually inject energy into the model and deteriorate performance
of the DAS over time. Similar problems may arise due to lack of correlation between some state
elements updated with the same transforms, so that even in presence of local observations the
ensemble spread for some elements may hardly reduce after assimilation, yet the ensemble anomalies
are inflated.
To avoid this behaviour EnKF-C currently restricts inflation by specified fraction (half by default)
of the the spread reduction factor calculated directly for each element of the state vector during
the update. For example, if inflation is specified as
INFLATION = 1.06 0.5
34

then the ensemble anomalies for any model state element will be inflated by 6 %, but no more than
1 + 0.5(σf /σa − 1), where σf and σa represent the forecast and analysis ensemble spreads for this
element. Specification
INFLATION = 1.06
is equivalent to
INFLATION = 1.06 1
Capping inflation by the magnitude of reduction of the ensemble spread is a default mode of
application of inflation in EnKF-C; to revert to the uniform inflation add qualifier PLAIN to the
entry INFLATION in the main parameter file, e.g.:
INFLATION = 1.06 PLAIN
The common inflation settings in the main parmater file can be overwritten by settings for particular
model variables specified in the model parameter file (sec. 2.4.2).

2.9

DA tuning

Following are the main parameters for DAS tuning in EnKF-C:
- R-factors;
- inflation magnitudes;
- localisation radii.
The R-factors can be defined for each observation type. They represent scaling coefficients for the
corresponding observation error variances and affect the impact of these observations: increasing
R-factor decreases the impact of observations and vice versa. Specifying R-factor equal k produces
the same increment as reducing the ensemble spread by k1/2 times.
The main parameter file defines the base R-factor common for all observation types. It is possible
to specify additional R-factors for observations of each type (sec. 2.4.4); the resulting R-factor for
an observation is then given by multiplication of the common R-factor and the additional R-factor
specified for observations of this type.
Multiplicative inflation can be seen as an additional forgetting factor in the KF. In EnKF-C one
can specify the inflation multiple for analysed ensemble anomalies, e.g.:
> grep INFLATION
INFLATION = 1.05
> grep temp -A 1
VAR = temp
INFLATION = 1.07

main.prm
PLAIN
model.prm
PLAIN

35

In this case all model variables except “temp” will have inflation of 5 %, while “temp” will have
inflation of 7 %. The ability to define different inflation rates for different variables can be useful
for non-dynamical variables, such as estimated biases, helping to avoid the ensemble collapse for
them. In general, to retain dynamical balances one should rather avoid using different inflation
magnitudes across model variables. Note that even small inflation can substantially affect the
ensemble spread established in the course of evolution of the system. By default EnKF-C applies
adaptive capping of inflation (sec. 2.8.1).
Localisation radius is defined by the entry LOCRAD in the parameter file. Specifically, this entry
defines the support localisation radius (in km). This is different to the “effective” localisation
radius, which is defined sometimes as e1/2 ≈ 1.65 - folding distance. For the Gaspary and Cohn’s
taper function used in EnKF-C the effective radius is approximately 3.5 times smaller than the
support radius.
Increasing the localisation radius increases the number of local observations and hence the overall
impact of observations. To compensate this in a system with horizontal localisation one has to
change the R-factor as the square of the localisation radius.

2.10

Point logs

It is often desirable to investigate sensitivity of model variables to observations or, more generally,
analyse certain features of the DAS and their behaviour over time. In practice such analysis can be
logistically complicated due to limitations on storage and/or access to it. Yet, it is usually feasible
to save this information in a number of specified locations.
EnKF-C provides capability of saving complete DA related information for specified horizontal
locations in so called “point logs”. The locations are specified in the main parameter file, e.g.:
POINTLOG
POINTLOG
POINTLOG
POINTLOG

94
78
57
86

134
111
51
191

Here the information will be saved for points with horizontal grid coordinates (94, 134), (78, 111),
(67, 51), and (86, 191) in files pointlog_94,134.nc, pointlog_78,111.nc, and so on. By default
the coordinates are defined on the first grid in the grid parameter file (sec. 2.4.3); otherwhile the
grid can be specified by entering the grid name as an additional parameter, .e.g:
POINTLOG 86 191 t-grid
Following is an example of the header of the saved point log file in NetCDF format:
netcdf pointlog_57\,51 {
dimensions:
m = 96 ;
p = 515 ;

36

nk-0 = 25 ;
nk-1 = 25 ;
nk-2 = 25 ;
variables:
int obs_ids(p) ;
float lcoeffs(p) ;
float lon(p) ;
float lat(p) ;
float depth(p) ;
float obs_val(p) ;
float obs_std(p) ;
float obs_fi(p) ;
float obs_fj(p) ;
float obs_fk(p) ;
int obs_type(p) ;
obs_type:SLA = 0 ;
obs_type:RFACTOR_SLA = 2. ;
obs_type:LOCRAD_SLA = 200. ;
obs_type:SST = 1 ;
obs_type:RFACTOR_SST = 4. ;
obs_type:LOCRAD_SST = 200. ;
obs_type:TEM = 2 ;
obs_type:RFACTOR_TEM = 8. ;
obs_type:LOCRAD_TEM = 800. ;
obs_type:SAL = 3 ;
obs_type:RFACTOR_SAL = 8. ;
obs_type:LOCRAD_TEM = 800. ;
float obs_date(p) ;
obs_date:units = "days from 5875.5 days since 1990-01-01" ;
float s(p) ;
float S(m, p) ;
double X5(m, m) ;
X5:long_name = "ensemble transform calculated for this (i,j) location in grid-0
(rho-grid)" ;
float X5_actual(m, m) ;
X5_actual:long_name = "the actual (interpolated) ensemble transform used at this
(i,j) location in grid-0 (rho-grid)" ;
float zeta(m) ;
float zeta_an(m) ;
zeta_an:INFLATION = 1.1f, 0.5f ;
float temp(nk-0, m) ;
float temp_an(nk-0, m) ;
temp_an:INFLATION = 1.1f, 0.5f ;
float salt(nk-0, m) ;
float salt_an(nk-0, m) ;
salt_an:INFLATION = 1.1f, 0.5f ;
float u(nk-1, m) ;
float u_an(nk-1, m) ;
u_an:INFLATION = 1.1f, 0.5f ;
float v(nk-2, m) ;
float v_an(nk-2, m) ;
v_an:INFLATION = 1.1f, 0.5f ;
float ubar(m) ;
float ubar_an(m) ;
ubar_an:INFLATION = 1.1f, 0.5f ;
float vbar(m) ;
float vbar_an(m) ;
vbar_an:INFLATION = 1.1f, 0.5f ;

37

// global attributes:
:date = "5875.5 days since 1990-01-01" ;
:i = 57 ;
:j = 51 ;
:lon = 151.75 ;
:lat = -39.849998 ;
:depth = 4632.623046875 ;

This data makes it possible to check DA algorithms by reproducing the ensemble transforms (for
EnKF) or weights (for EnOI) calculated by EnKF-C from S and s according to section 1.3.3; restore
observations from s by using the corresponding R-factors and localisation coefficients; to monitor
the ensemble spread for each model variable; calculate inflation applied to the analysed anomalies;
calculate impacts of particular observations; and so on.

2.11

Use of innovation statistics for model validation

EnKF-C can calculate innovation statistics for validating a model against observations only, without
data assimilation. The pre-requisites are (i) observations and (ii) model dump readable by the code,
and possibly (iii) auxiliary files for projecting the model state to observation space (e.g. grid specs
and mean SSH). To get the innovation statistics one needs to:
- set up the parameter files in a normal way (MODE = EnOI), omitting the ensemble directory
and assimilation related parameters;
- run enkf_prep;
- run enkf_calc with additional parameter --forecast-stats-only.
The results will be written to the log of enkf_calc. An example of using this functionality is
available by running make stats in examples/1 (see sec. 2.3).

2.12

Bias correction

It is possible to estimate and correct bias for a model variable with the EnKF by generating and
using an ensemble of bias fields. These bias fields need to be subtracted from the corresponding
observation forecasts. This is accomplished by specifying a secondary variable in the observation
type entry VAR and by passing the name of this variable to the corresponding observation functions,
which need to take care for subtracting the bias from the model forecasts.
Because bias fields are usually assumed to persist (not change) during propagation, one may need
to make specific settings for inflation for them to avoid their collapse (loss of spread) over time.
Another possibility is to introduce a “forgetting” stochastic model for bias fields, for example:
xi+1 = λ xi + (1 − λ2 )1/2 σ,
where α is the forgetting factor, 0 < λ < 1, 1 − λ ≪ 1, and σ ∼ σ0 N (0, 1), where σ0 is the standard
deviation of the error of x. This can be specified for a model variable via entry RANDOMISE in the
corresponding section of the model parameter file (sec. 2.4.2).
38

2.13

Grids

Currently EnKF-C makes three main general assumptions about grids:
- the vertical grid coordinates correspond to layer centres, either for Z or sigma layers;
- the horizontal coordinates are geographic (longitude, latitude);
- the order of vertical coordinate (top to bottom or bottom to top) corresponds to that of 3D
model fields.
If the available model grids do not meet the above assumptions, one needs to generate the corresponding auxiliary grids and use them instead.
Further, the horizontal grid must be either rectangular in longitude-latitude, with constant or
irregular step; or represent a structured simply connected quadrilateral grid, referred to elsewhere
as “curvilinear”. In the first case mapping from physical to grid space is done independently by
X and Y coordinates; in the second case the grid is rendered into spatial binary tree to achieve
logarithmic search time. Note that some curvilinear grids used in practice (e.g. tri-polar ORCA
grids) are not simply connected in grid space. Such grids can be used in EnKF-C as described in
section 2.4.3.

2.14

System issues

2.14.1

Memory footprint

To reduce the memory usage, most of the potentially big arrays in EnKF-C use float data type.
The memory footprint of prep is defined by the size of the measurement structure and the number
of observations. It is rarely a problem.
For calc it is mainly defined by the size of ensemble observation anomalies, which require p × m × 4
bytes for storage. For example, with 3 · 106 superobservations and 102 ensemble members the size
of this array would be about 1.2 GB per CPU, which should be manageable on most contemporary
systems. If the footprint becomes too big, one may consider reducing the number of observations
by a coarser superobing (setting the parameter SOBSTRIDE to 2 or more) or reducing the number
of cores per node used.
The footprint of update is mainly defined by the size of the array of simultaneously updated horizontal fields. For example, for a 1500 × 3600 horizontal field, 100 ensemble members and simultaneous
update of 2 fields the size of this array would be about 4.3 GB per CPU. It can be reduced to
2.15 GB by setting the parameter FIELDBUFFERSIZE to 1. Note that reducing the number of simultaneously updated fields specified by parameter FIELDBUFFERSIZE increases I/O (the X5 array
is read from disk Nf /Nb times, where Nf is the total number of horizontal fields, and Nb is the
number of simultaneously updated fields) and reduces computational effectiveness (the X5 array
needs to be interpolated horizontally Nf /Nb times).

39

2.14.2

Exit action

When exiting on an error, EnKF-C by default prints the stack trace, which allows to trace the
exit location in the code. Another option – to generate a segmentation fault – can be activated by
setting EXITACTION = SEGFAULT in the parameter file. Note that when run on multiple processors,
this can result in segmentation faults on more than one processor (but not necessarily on every
engaged processor, as some processes can also be forced to exit by MPI_abort()). If the system is
set to generate core dumps, they can indeed be used for investigating the final state of the program.

2.14.3

Dependencies and compilation issues

Compiling EnKF-C requires the following external libraries:
- netcdf;
- lapack (or mkl rt);
- openmpi;
- gridutils.
EnKF-C also relies on qsort_r(), which may be lacking in older systems. In such cases use compile
flag -DINTERNAL_QSORT_R to activate the internal version of this procedure.
Notes:
1. Using Intel’s version of Lapack library – Intel Math Kernel Library – can improve performance
over Lapack compiled with gfortran.
2. Grid utilities (libgu) is necessary for handling curvilinear grids only.
-DNO_GRIDUTILS to compile without it.

Add compile flag

2.15

Possible problems

2.15.1

calc becomes too slow after increasing localisation radius

This is due to the increased number of local observations.
The local observations are sought by using so called k-d tree. Normally it works well, but can
become a bottleneck when the number of local observations becomes very large. If this happens,
one may use a coarser superobing to reduce the number of assimilated observations, by increasing
parameter SOBSTRIDE from the default value of 1 to 2 or more.

40

Acknowledgements
EnKF-C has been developed during author’s work with Bureau of Meteorology on Bluelink project.
The author used his knowledge of TOPAZ (Sakov et al., 2012) and BODAS (Oke et al., 2008) codes
and borrowed from them a number of design solutions and features. Paul Sandery was the first
user of this code (apart from the author), and his enthusiastic support is cheerfully acknowledged.

41

References
Andrews, A., 1968: A square root formulation of the Kalman covariance equations. AIAA J., 6,
1165–1168.
Bishop, C. H., B. Etherton, and S. J. Majumdar, 2001: Adaptive sampling with the ensemble
transform Kalman filter. part I: theoretical aspects. Mon. Wea. Rev., 129, 420–436.
Evensen, G., 1994: Sequential data assimilation with a nonlinear quasi-geostrophic model using
Monte-Carlo methods to forecast error statistics. J. Geophys. Res., 99, 10143–10162.
— 2003: The Ensemble Kalman Filter: theoretical formulation and practical implementation.
Ocean Dynam., 53, 343–367.
— 2004: Sampling strategies and square root analysis schemes for the EnKF. Ocean Dynam., 54,
539–560.
Evensen, G. and P. J. van Leeuwen, 2000: An ensemble Kalman smoother for nonlinear dynamics.
Mon. Wea. Rev., 128, 1852–1867.
Gaspari, G. and S. E. Cohn, 1999: Construction of correlation functions in two and three dimensions. Q. J. R. Meteorol. Soc., 125, 723–757.
Hamill, T. M. and J. S. Whitaker, 2001: Distance-dependent filtering of background error covariance
estimates in an ensemble Kalman filter. Mon. Wea. Rev., 129, 2776–2790.
Houtekamer, P. L. and H. L. Mitchell, 2001: A sequential ensemble Kalman filter for atmospheric
data assimilation. Mon. Wea. Rev., 129, 123–137.
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 filtering. Tellus, 56A, 273–277.
Hunt, B. R., E. J. Kostelich, and I. Szunyogh, 2007: Efficient data assimilation for spatiotemporal
chaos: A local ensemble transform Kalman filter. Physica D , 230, 112–126.
Kalman, R. E., 1960: A new approach to linear filtering and prediction problems. J. Basic. Eng.,
82, 35–45.
Oke, P. R., G. B. Brassington, D. A. Griffin, and A. Schiller, 2008: The Bluelink ocean data
assimilation system (BODAS). Ocean Model., 21, 46–70.
Ott, E., B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E. Kalnay, D. J.
Patil, and J. A. Yorke, 2003, rev. 2005: A local ensemble Kalman filter for atmospheric data
assimilation. http://arxiv.org/abs/physics/0203058 .
42

Sakov, P. and L. Bertino, 2011: Relation between two common localisation methods for the EnKF.
Comput. Geosci., 15, 225–237.
Sakov, P., F. Counillon, L. Bertino, K. A. Lisæter, P. R. Oke, and A. Korablev, 2012: TOPAZ4:
an ocean-sea ice data assimilation system for the North Atlantic and Arctic. Ocean Science, 8,
633–656.
Sakov, P., G. Evensen, and L. Bertino, 2010: Asynchronous data assimilation with the EnKF.
Tellus, 62A, 24–29.
Sakov, P. and P. R. Oke, 2008a: A deterministic formulation of the ensemble Kalman filter: an
alternative to ensemble square root filters. Tellus, 60A, 361–371.
— 2008b: Implications of the form of the ensemble transformations in the ensemble square root
filters. Mon. Wea. Rev., 136, 1042–1053.
Verlaan, M. and A. W. Heemink, 1997: Tidal flow forecasting using reduced rank square root filters.
Stoch. Hydrol. Hydraul., 11, 349–368.

43

Abbreviations
CL
DEnKF
DA
DAS
DAW
DFS
EKF
EnKF
EnOI
ETKF
ETM
FGAT
KF
KS
LA
SDAS
SRF
SVD

-

covariance localisation
deterministic EnKF
data assimilation
data assimilation system
data assimilation window
degrees of freedom of signal
extended Kalman filter
ensemble Kalman filter
ensemble optimal interpolation
ensemble transform Kalman filter
ensemble transform matrix
first guess at appropriate time
Kalman filter
Kalman smoother
local analysis
state of data assimilation system
spread reduction factor
singular value decomposition

44

Symbols
General symbols
x (small, bold)
1
0
A (capital, bold)
I
U
AT
A1/2
tr(A)
H ◦ M(x)
A◦B

-

a vector
a vector with all elements equal to 1
a vector with all elements equal to 0
a matrix
an identity matrix
a unitary matrix, UUT = I
transposed matrix A
the unique positive definite square root of a positive definite matrix A
trace of A
H [M(x)]
by-element, or Hadamard, or Schur product of matrices

45

DA related symbols
m
n
p
A
E
G
H
H
J
M
M
P
Q
R
S
TL
TR
Up
X5
s
x
y
w
(·)f
(·)a
(·)i

-

ensemble size
state size
number of observations
ensemble anomalies, A = E − x1T
ensemble
an intermediate matrix in the EnKF analysis, G ≡ (I + ST S)−1 ST = ST (I + SST )−1
nonlinear observation operator; in linear case – affine observation operator
linearised observation operator, H = ∇H(x)
cost function
nonlinear model operator; in linear case – affine model operator
linearised model operator, M = ∇M(x)
state error covariance estimate; also used as abbreviation for AAT /(m − 1)
model error covariance
observation error covariance
√
normalised ensemble observation anomalies, S = R−1/2 HA/ m − 1
left-multiplied ensemble transform matrix, Aa = TL Af
right-multiplied ensemble transform matrix, Aa = Af TR
a unitary mean-preserving matrix, Up (Up )T = I, Up 1 = 1
historic symbol for the full ensemble
transform
matrix,
E a = E f X5


√
normalised innovation, s = R−1/2 y − H(xf ) / m − 1
state estimate
observation vector
vector of linear coefficients for updating the mean, xa = xf + Af w
forecast expression
analysis expression
either expression at cycle i or ith element of a vector

(·)

-

local expression for state element i

-

local expression for observation o

i

{o}

(·)

46

