Koren B., Vuik K. (Editors) Advanced Computational Methods in Science and Engineering
Подождите немного. Документ загружается.
Parallel Scientific Computing on Loosely Coupled Networks of Computers
Fig. 8 Artist’s impression of fish with Immersed Boundaries (artwork by Tobias Baanders).
References
1. David P. Anderson, Jeff Cobb, Eric Korpela, Matt Lebofsky, and Dan Werthimer.
SETI@home: an experiment in public–resource computing. Commun. ACM, 45(11):56–61,
2002.
2. Owe Axelsson. Iterative Solution Methods. Cambridge University Press, New York, NY,
USA, 1994.
3. Jacques M. Bahi, Sylvain Contassot-Vivier, and Rapha¨el Couturier. Evaluation of the asyn-
chronous iterative algorithms in the context of distant heterogeneous clusters. Parallel Com-
put., 31(5):439–461, 2005.
4. D. El Baz. A method of terminating asynchronous iterative algorithms on message passing
systems. Parallel Algorithms and Applications, 9:153–158, 1996.
5. Didier El Baz, Pierre Spiteri, Jean Claude Miellou, and Didier Gazen. Asynchronous iterative
algorithms with flexible communication for nonlinear network flow problems. J. Parallel
Distrib. Comput., 38(1):1–15, 1996.
6. Micah Beck, Dorian Arnold, Alessandro Bassi, Fran Berman, Henri Casanova, Jack Dongarra,
Terry Moore, Graziano Obertelli, James Plank, Martin Swany, Sathish Vadhiyar, and Rich
Wolski. Middleware for the use of storage in communication. Parallel Comput., 28(12):1773–
1787, 2002.
103
Tijmen P. Collignon and Martin B. van Gijzen
7. Michele Benzi. Preconditioning techniques for large linear systems: a survey. J. Comput.
Phys., 182(2):418–477, 2002.
8. Fran Berman, Geoffrey Fox, and Anthony J. G. Hey. Grid Computing: Making the Global
Infrastructure a Reality. John Wiley & Sons, Inc., New York, NY, USA, 2003.
9. Dimitri P. Bertsekas and John N. Tsitsiklis. Parallel and Distributed Computation: Numerical
Methods. Prentice–Hall, Englewood Cliffs, N.J., 1989. republished by Athena Scientific,
1997.
10. Rob H. Bisseling. Parallel Scientific Computation: A Structured Approach Using BSP and
MPI. Oxford University Press, 2004.
11.
˚
A. Bj¨orck. Solving linear least squares problems by Gram–Schmidt orthogonalization. BIT,
7:1–21, 1967.
12. Kostas Blathras, Daniel B. Szyld, and Yuan Shi. Timing models and local stopping criteria for
asynchronous iterative algorithms. Journal of Parallel and Distributed Computing, 58(3):446–
465, 1999.
13. T. Brady, E. Konstantinov, and A. Lastovetsky. SmartNetSolve: High level programming
system for high performance Grid computing. Rhodes Island, Greece, 25-29 April 2006 2006.
IEEE Computer Society. CD-ROM/Abstracts Proceedings.
14. Eddy Caron, Bruno Del-Fabbro, Fr´ed´eric Desprez, Emmanuel Jeannot, and Jean-Marc Nicod.
Managing data persistence in network enabled servers. Sci. Program., 13(4):333–354, 2005.
15. Eddy Caron and Fr´ed´eric Desprez. DIET: A scalable toolbox to build network enabled servers
on the Grid. International Journal of High Performance Computing Applications, 20(3):335–
352, 2006.
16. A. T. Chronopoulos and C. W. Gear. S–step iterative methods for symmetric linear systems.
J. Comput. Appl. Math., 25(2):153–168, 1989.
17. Tijmen P. Collignon and Martin B. van Gijzen. Implementing the Conjugate Gradient Method
on a grid computer. In Proceedings of the International Multiconference on Computer Science
and Information Technology, Volume 2, October 15–17, 2007, Wisla, Poland, pages 527–540,
2007.
18. Tijmen P. Collignon and Martin B. van Gijzen. Solving large sparse linear systems efficiently
on Grid computers using an asynchronous iterative method as a preconditioner. Technical
report, Delft University of Technology, Delft, the Netherlands, 2008. DUT report 08–08.
19. Tijmen P. Collignon and Martin B. van Gijzen. Two implementations of the preconditioned
Conjugate Gradient method on a heterogeneous computing grid with applications to 3D bub-
bly flow problems. Technical report, Delft University of Technology, Delft, the Netherlands,
2008. DUT report 08–xx.
20. Rapha¨el Couturier, Christophe Denis, and Fabienne J´ez´equel. GREMLINS: a large sparse
linear solver for grid environment. Parallel Computing, December 2008.
21. Rapha¨el Couturier and St´ephane Domas. CRAC: a Grid Environment to solve Scientific Ap-
plications with Asynchronous Iterative Algorithms. In 21th IEEE and ACM Int. Symposium
on Parallel and Distributed Processing Symposium, IPDPS’2007, page 289 (8 pages), Long
Beach, USA, March 2007. IEEE computer society press.
22. J. Daniel, W. B. Gragg, L. Kaufman, and G. W. Stewart. Reorthogonalization and stable
algorithms for updating the Gram–Schmidt QR factorization. Mathematics of Computation,
30:772–795, 1976.
23. Eric de Sturler. Truncation strategies for optimal Krylov subspace methods. SIAM J. Numer.
Anal., 36(3):864–889, 1999.
24. Frederic Desprez and Emmanuel Jeannot. Improving the GridRPC model with data persis-
tence and redistribution. In ISPDC ’04: Proceedings of the Third International Symposium
on Parallel and Distributed Computing/Third International Workshop on Algorithms, Models
and Tools for Parallel Computing on Heterogeneous Networks (ISPDC/HeteroPar’04), pages
193–200, Washington, DC, USA, 2004. IEEE Computer Society.
25. K.D. Devine, E.G. Boman, R.T. Heaphy, R.H. Bisseling, and U.V. Catalyurek. Parallel hy-
pergraph partitioning for scientific computing. In Proc. of 20th International Parallel and
Distributed Processing Symposium (IPDPS’06). IEEE, 2006.
104
Parallel Scientific Computing on Loosely Coupled Networks of Computers
26. J. Dongarra and A. Lastovetsky. An overview of heterogeneous high performance and Grid
computing. Engineering the Grid: Status and Perspective, February 2006 2006.
27. Jack Dongarra, Ian Foster, Geoffrey Fox, William Gropp, Ken Kennedy, Linda Torczon, and
Andy White, editors. Sourcebook of Parallel Computing. Morgan Kaufmann, 2003.
28. Jack Dongarra, Yinan Li, Zhiao Shi, Don Fike, Keith Seymour, and Asim YarKhan. Homepage
of NetSolve/GridSolve, 2007. .
29. Stanley C. Eisenstat, Howard C. Elman, and Martin H. Schultz. Variational iterative methods
for nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20:345–357, 1983.
30. Folding. Folding@home distributed computing. .
31. I. Foster and C. Kesselman. The Grid: Blueprint for a new Computing Infrastructure. Morgan
Kaufman Publishers, second edition, 2004.
32. J. Frank and C. Vuik. On the construction of deflation–based preconditioners. SIAM J. Sci.
Comput., 23(2):442–462, 2001.
33. Andreas Frommer and Daniel B. Szyld. Asynchronous iterations with flexible communication
for linear systems. Calculateurs Parall`eles R´eseaux et Syst`emes R´epartis, 10:421–429, 1998.
34. Gene H. Golub and Charles F. Van Loan. Matrix Computations (Johns Hopkins Studies in
Mathematical Sciences). The Johns Hopkins University Press, October 1996.
35. Magnus R. Hestenes and Eduard Stiefel. Methods of Conjugate Gradients for solving linear
systems. Journal of Research of National Bureau Standards, 49:409–436, 1952.
36. A. Lastovetsky, X. Zuo, and P. Zhao. A non–intrusive and incremental approach to enabling
direct communications in RPC–based grid programming systems. Technical report, 2006.
37. Craig Lee, Hidemoto Nakada, and Yusuke Tanimura. GridRPC Working Group, 2007.
.
38. J. C. Miellou, D. El Baz, and P. Spiteri. A new class of asynchronous iterative algorithms with
order intervals. Math. Comput., 67(221):237–255, 1998.
39. Y. Notay. Flexible Conjugate Gradients. SIAM Journal on Scientific Computing, 22:1444–
1460, 2000.
40. Y. Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied
Mathematics, Philadelphia, PA, USA, 2003.
41. Youcef Saad. A flexible inner–outer preconditioned GMRES algorithm. SIAM J. Sci. Comput.,
14(2):461–469, 1993.
42. Youcef Saad and Martin H. Schultz. GMRES: a generalized minimal residual algorithm for
solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput., 7(3):856–869, 1986.
43. Mitsuhisa Sato, Taisuke Boku, and Daisuke Takahashi. OmniRPC: a Grid RPC system for par-
allel programming in cluster and Grid environment. In CCGRID ’03: Proceedings of the 3rd
International Symposium on Cluster Computing and the Grid, pages 206–213, Washington,
DC, USA, 2003. IEEE Computer Society.
44. K Seymour, A YarKhan, S Agrawal, and J Dongarra. NetSolve: Grid enabling scientific com-
puting environments. In L. Grandinetti, editor, Grid Computing and New Frontiers of High
Performance Processing. Elsevier, 2005.
45. Keith Seymour, Hidemoto Nakada, Satoshi Matsuoka, Jack Dongarra, Craig Lee, and Henri
Casanova. Overview of GridRPC: A Remote Procedure Call API for Grid Computing. In
GRID ’02: Proceedings of the Third International Workshop on Grid Computing, pages 274–
278, London, UK, 2002. Springer–Verlag.
46. Valeria Simoncini and Daniel B. Szyld. Flexible inner–outer Krylov subspace methods. SIAM
J. Numer. Anal., 40(6):2219–2239, 2002.
47. Valeria Simoncini and Daniel B. Szyld. Recent computational developments in Krylov sub-
space methods for linear systems. Numerical Linear Algebra with Applications, 14:1–59,
2007.
48. Barry F. Smith, Petter E. Bjørstad, and William Gropp. Domain Decomposition: Parallel
Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press,
Cambridge, 1996.
49. Peter Sonneveld and Martin B. van Gijzen. IDR(s): a family of simple and fast algorithms for
solving large nonsymmetric linear systems. Technical report, Delft University of Technology,
Delft, the Netherlands, 2007. DUT report 07–07.
105
Tijmen P. Collignon and Martin B. van Gijzen
50. StarPlane. Application-specific management of photonic networks, 2007.
.
51. Thomas Sterling, Ewing Lusk, and William Gropp, editors. Beowulf Cluster Computing with
Linux. MIT Press, Cambridge, MA, USA, 2003.
52. Y. Tanaka, H. Nakada, S. Sekiguchi, T. Suzumura, and S. Matsuoka. Ninf–G: A reference
implementation of RPC–based programming middleware for Grid computing. Journal of
Grid Computing, 1(1):41–51, 2003.
53. James D. Teresco, Karen D. Devine, and Joseph E. Flaherty. Numerical Solution of Partial
Differential Equations on Parallel Computers, chapter Partitioning and Dynamic Load Bal-
ancing for the Numerical Solution of Partial Differential Equations. Springer–Verlag, 2005.
54. A. Toselli and O. B. Widlund. Domain Decomposition: Algorithms and Theory, volume 34.
Springer Series in Computational Mathematics, Springer, Berlin, Heidelberg, 2005.
55. H. A. van der Vorst. Bi–CGSTAB: A fast and smoothly converging variant of Bi–CG for the
solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Comput-
ing, 13(2):631–644, 1992.
56. H.A. van der Vorst and C. Vuik. GMRESR: a family of nested GMRES methods. Num. Lin.
Alg. Appl., 1(4):369–386, 1994.
57. Henk A. van der Vorst. Iterative Krylov Methods for Large Linear systems. Cambridge Uni-
versity Press, Cambridge, 2003.
58. Brendan Vastenhouw and Rob H. Bisseling. A two–dimensional data distribution method for
parallel sparse matrix-vector multiplication. SIAM Rev., 47(1):67–95, 2005.
59. Flavien Vernier, Jacques M. Bahi, Sylvain Contassot-Vivier, and Raphael Couturier. A decen-
tralized convergence detection algorithm for asynchronous parallel iterative algorithms. IEEE
Trans. Parallel Distrib. Syst., 16(1):4–13, 2005.
60. P. Wesseling. An Introduction to Multigrid Methods. John Wiley & Sons, Chichester, 1992.
61. Samuel Williams, John Shalf, Leonid Oliker, Shoaib Kamil, Parry Husbands, and Katherine
Yelick. The potential of the cell processor for scientific computing. In CF ’06: Proceedings
of the 3rd conference on Computing frontiers, pages 9–20, New York, NY, USA, 2006. ACM.
62. Asim YarKhan, Keith Seymour, Kiran Sagi, Zhiao Shi, and Jack Dongarra. Recent devel-
opments in GridSolve. International Journal of High Performance Computing Applications
(IJHPCA), 20(1):131–141, 2006.
63. X. Zuo and A. Lastovetsky. Experiments with a software component enabling NetSolve with
direct communications in a non–intrusive and incremental way. In Proceedings of the 21st
International Parallel and Distributed Processing Symposium (IPDPS 2007), Long Beach,
California, USA, 26-30 March 2007 2007. IEEE Computer Society.
106
Data Assimilation Algorithms for Numerical
Models
A.W. Heemink, R.G. Hanea, J. Sumihar, M. Roest, N. Velzen and M. Verlaan
Abstract To understand and predict the behavior of a system one can use measure-
ments or one can develop physically based numerical models. In many applications
however neither of these approaches is able to provide an accurate description of
the dynamic behavior of the system. A model is always a simplification of the real
world while measurements seldom produce a complete picture of the system be-
haviour. Using data assimilation techniques measurements and model results are
both used to obtain an optimal estimate of the state of the system. In this article we
present an overview of methods available to assimilate data into a numerical model.
Attention is concentrated on variational methods and on Kalman filtering. The main
problem of using these advanced data assimilation schemes is the huge computa-
tional burden that is required for solving real life problems. For variational methods
the adjoint model implementation is essential to obtain an efficient data assimilation
algorithm. For Kalman filtering problems a number of approximate algorithms have
been introduced recently: Ensemble Kalman filters and Reduced Rank filters. These
algorithms make the application of Kalman filtering to large-scale data assimilation
problems feasible. After a brief introduction to the most important data assimila-
tion approaches we will discuss the advantages and disadvantages of the various
methods and present a number of real life applications.
1 Introduction
Measurements can be used to develop statistical models for predicting the behavior
of environmental processes. However these types of models are derived from the
data and do not include physical knowledge of the process. Furthermore, measure-
ments alone do generally not provide a complete picture of the process. Especially
in case of processes that vary in space and time it is very hard to reconstruct the
A.W. Heemink
Delft University of Technology, e-mail: A.W.Heemink@tudelft.nl
107
Lecture Notes in Computational Science and Engineering 71,
DOI 10.1007/978-3-642-03344-5_5,
B. Koren and C. Vuik (eds.), Advanced Computational Methods in Science and Engineering,
© Springer-Verlag Berlin Heidelberg 2010
A.W. Heemink et al.
spatial and temporal patterns only from data. Physically based numerical models
produce results that are spatially and temporally consistent. However these models
are usually not able to accurately reproduce the measurements that are available.
The information provided by the models and by the measurement information is
often complementary. Therefore it is important to study a methodology for integrat-
ing measurements and physically based mathematical models. This methodology is
called data assimilation. By using numerical models that are based on physical laws
and that are continuously adapted by the measurements available the two sources of
information of the process, model information and measurement information, can
be integrated.
Data assimilation can be defined as a procedure to incorporate data into a model
simulation so as to improve the predictions. However, assimilating data into a nu-
merical model is far from trivial. The simplest data assimilation procedure is to
overwrite the model values at the measurement locations with the observed data.
Inserting the data in this way into a numerical model is in general not a satisfactory
method. It leaves the model dynamically unbalanced and introduces spurious waves
into the model. These short waves may even cause instabilities of the underlying
numerical model.
A common data assimilation technique used in numerical weather prediction is
optimal interpolation. Here some estimates of the error statistics of the numerical
model are used to correct the results of the model using the measurements. How-
ever, since these error statistics have to be determined by adopting some ad-hoc
statistical assumptions, the correction produced by optimal interpolation is again
not consistent with the underlying numerical model. As a consequence the use of
optimal interpolation often still yields unrealistic correction or instabilities.
More accurate data assimilation methods are variational data assimilation and
Kalman filtering. The basic idea of these data assimilation methods is to use the data
to only correct the weak points in the model. Weak parts of the model may be due to
uncertainty in initial and boundary conditions or imperfectly known model param-
eters. The data is not allowed to modify the accurate parts of the model. As a result
these types of data assimilation problems are in fact, inverse problems. The specified
inputs (model uncertainties) have to be reconstructed from the output (measure-
ments). A variational approach or Kalman filtering solves these inverse problems
accurately. For linear problems it can be shown that both approaches produce ex-
actly the same results for the same problem formulation. Optimal interpolation does
not solve an inverse problem. It produces corrections for the model output without
reconstructing model uncertainties. As a result the variational method and Kalman
filtering are superior to optimal interpolation.
In the last decennium the variational approach and Kalman filtering have gained
acceptance as powerful frameworks for data assimilation. However, both methods
require a very large computational burden, at least an order of magnitude larger than
the computational effort required for the underlying numerical model. This is the
main disadvantage of these methods compared to optimal interpolation that requires
only a small increase in computer time.
108
Data Assimilation Algorithms for Numerical Models
Starting point for the data assimilation methodology is a state space representa-
tion of the numerical model and the measurements. Let us assume that modeling
techniques have provided us with a deterministic state space representation of the
form:
X
k+1
= f (X
k
,k) + B(k)u
k
, X
0
= x
0
. (1)
Here the X
k
is the system state, u
k
is the input of the system, f is a nonlinear function,
and B(k) is an input matrix. For a numerical model that describes the behavior of a
physical process in space and time, the state consists of all the variables in all the
grid points of the model at a certain time, while the function f in this case represents
one time step of the numerical scheme of the model.
The measurements taken from the actual system are assumed to be available
according to the relation:
Z
k
= m(X
k
,k), (2)
where Z
k
is a vector containing the measurements and m is a nonlinear function that
specifies the relation between the model results and the measurements.
In this chapter we describe in Section 2 the basic idea of the variational approach
and discuss a number of extensions. In Section 3 we introduce the Kalman filter as
data assimilation framework. Here we present a number of filter algorithms for solv-
ing large-scale data assimilation problems. In Section 4 a software environment for
data assimilation is presented and the advantages of having a general framework for
applying different algorithms for different applications. Some of the large scale real
life applications are presented in the following chapters. In Section 5 applications
in coastal sea modeling are presented, in the specific cases of storm surge predic-
tion and assimilation of high frequency (HF) radar data into a coastal ocean model.
A large scale atmospheric chemistry application is presented in Section 6 together
with a comparison of the performance of different sequential algorithms.
2 Variational data assimilation
2.1 Data assimilation formulated as a minimization problem
If it can be assumed that the only uncertainties of the model (1)-(2) are introduced by
a number of poorly known parameters, the data assimilation problem can be formu-
lated as a deterministic parameter estimation problem. Rewrite the model according
to:
X
k+1
= f (X
k
, p, k ) + B(k, p)u
k
, X
0
= x
0
, (3)
Z
k
= m(X
k
,k),
where p is a vector containing the uncertain parameters. Uncertain parameters may
be model parameters, initial conditions or inputs. The state X
k
is now dependent on
the parameters X
k
(p). In order to estimate the parameters we first define a criterion
109
A.W. Heemink et al.
J(p) as a measure for the distance between the measurements and the model results:
J(p) =
K
∑
k=1
Z
k
−m(X
k
(p),k)
T
R(k)
−1
Z
k
−m(X
k
(p),k)
. (4)
Here the generalized least squares criterion has been chosen to define J. The covari-
ance matrix R(k) is a weighting matrix that takes into account the errors associated
with the measurements. This formulation is used very often in practice. The optimal
parameter p is found by minimizing the criterion J(p).
Prior information about the parameter values p
0
can be included by adding a
regularization term to the criterion:
J(p) =
K
∑
k=1
Z
k
−m(X
k
(p),k)
T
R(k)
−1
Z
k
−m(X
k
(p),k)
+(p − p
0
)
T
P
−1
0
(p − p
0
). (5)
Here P
0
is the covariance matrix of the prior information p
0
, modeling the uncer-
tainty associated with this prior information. Usually p
0
is the first guess of the
uncertain parameters and the starting value for the optimization procedure. The reg-
ularization term in the criterion prevents that parameter estimates become unrealis-
tic if the measurement information is limited. In this case the estimates will simply
remain close to the first guess (as they should be).
2.2 The adjoint model
In most practical data assimilation problems the number of uncertain parameters is
very large. For example, if p represents the initial condition of a numerical model,
the number of parameters is equal to the dimension of the system state. Most op-
timization techniques are not able to solve these very large-scale problems. In this
case the only feasible approach is to use a gradient-based optimization method and
to use a variational approach to efficiently compute the gradient of the criterion.
For every parameter p the system state X
k
has to satisfy the model constraint (3).
Therefore we can rewrite the criterion (4) as:
J(p) =
K
∑
k=1
Z
k
−m(X
k
(p),k)
T
R(k)
−1
Z
k
−m(X
k
(p),k)
+
K−1
∑
k=0
L
T
k+1
X
k+1
− f (X
k
, p, k ) −B(k, p)u
k
, (6)
where L
k
represents a Lagrange multiplier or adjoint state. Note that expression (6)
holds for any choice of L
k
. It is easy to show that if L
k
satisfies the system of adjoint
equations:
110
Data Assimilation Algorithms for Numerical Models
L
T
k
= L
T
k+1
∂ f (X
k
, p, k )
∂ X
k
+ 2
Z
k
−m(X
k
,k)
T
R(k)
−1
∂ m(X
k
,k)
∂ X
k
, (7)
for k = K −1, K −2, . .. , 1, with end condition L
K
= 0, the gradient of the criterion
can be computed by using the very simple expression:
∂ J
∂ p
= −
K−i
∑
k=0
L
T
k+1
∂ f (X
k
, p, k )
∂ p
+
∂ B(k, p)
∂ p
u
k
. (8)
Using this expression the gradient of the criterion can be determined from the results
of one forward simulation to evaluate the terms:
∂ f (X
k
, p, k )
∂ p
+
∂ B(k, p)
∂ p
u
k
, (9)
for k = 1, 2, ... , K, and the results L
k
, for k = K, K −1, ..., 1 of one adjoint
model simulation. The computational effort required to determine the gradient is
almost independent of the number of parameters. Therefore, by combining this idea
with a gradient based optimization scheme, a very efficient implementation can be
obtained, especially for data assimilation problems with many uncertain parameters.
This approach is often called variational data assimilation or the adjoint method.
2.3 Discussion
Variational data assimilation schemes are iterative schemes to minimize a criterion
J(p). The number of iterations needed depends on the optimization scheme used.
In practice quasi-Newton schemes are often used for data assimilation problems.
For these types of schemes it is known that if the criterion is exactly quadratic, the
iteration process converges in d
p
+ 1 steps, with d
p
as the dimension of p. In most
practical data assimilation problems the actual number of iterations can be chosen
significantly less since the largest improvement is often obtained in the first few
iterations.
One problem with the application of variational data assimilation schemes is that
the results of one complete forward simulation have to be stored. This is for most
real life data assimilation applications a serious problem. Therefore one generally
only stores the system state at a limited number of time steps. If results at interme-
diate time steps are necessary the states are recomputed using the forward model
again. This reduces the storage problem at the expense of additional computations.
Another problem with variational data assimilation is the implementation of the
adjoint model. Although the mathematical derivation of this model from the origi-
nal forward model is straightforward, the coding of the adjoint model can be very
difficult. Especially in case when the forward model has been developed step by
step over a long period by different programmers and many parts of the code are not
well understood (although they produce reliable results). Moreover, forward models
111
A.W. Heemink et al.
are often still under development and are being improved continuously. As a result
the adjoint should also be updated continuously. Recently adjoint compilers have
become available that produce adjoint code from the code of the forward model
(Kaminski et al., 2003). However, for complicated forward models these compil-
ers still produce unsatisfactory results. Probably the best way to develop an adjoint
implementation is to redesign and rewrite the forward code such that the adjoint
compiler is able to generate the adjoint. An advantage of this approach is that modi-
fications in the forward model can easily be included in the adjoint implementation.
Another approach to variational data assimilation with a comparable computa-
tional efficiency is based on model reduction (Ravindran, 2002, Vermeulen et al.,
2004, 2006). This approach does not require the implementation of the adjoint of the
tangent linear approximation of the original model. Using an ensemble of forward
model simulations an approximation of the covariance matrix of the model variabil-
ity is determined. A limited number of leading eigenvectors (EOF’s) of this matrix
are selected to define a model subspace. By projecting the original model onto this
subspace an approximate linear model is obtained. Once this reduced model is avail-
able, its adjoint can be implemented very easily and the minimization process can be
solved completely in reduced space with negligible computational costs. If neces-
sary, the procedure can be repeated a few times by generating new ensembles close
to the most recent estimate of the parameters.
The approach described in this section is a strong constraint variational method.
This means that the only uncertainties in the model are the specified parameters. The
model is considered to be perfect. A weak constraint variational approach allows for
model uncertainties too (Bennett, 2002, Heemink and Metzelaar, 1995 ). Like with
Kalman filtering model uncertainties are modeled by including stochastic forcing
term in the model equations. Algorithms based on a weak constraint formulation
are (even) more time consuming then the strong constraint approach and, therefore,
are not used very often yet.
3 Kalman filtering
3.1 The linear Kalman filter
In Section 2 model uncertainties were only due to poorly known parameters. Uncer-
tainties in a mathematical model can also be modeled by embedding the model in
a stochastic environment. For a linear system a stochastic state space representation
can be formulated as follows:
X
k+1
= F(k)X
k
+ B(k)u
k
+ G(k)W
k
, (10)
Z
k
= M(k)X
k
+V
k
. (11)
A Gaussian system noise process W
k
with zero mean and covariance matrix Q(k) is
introduced to take into account model uncertainties. G(k) is the noise input matrix.
112