Error subspace filtering for atmospheric chemistry data assimilation modeling

(1)

(2)

(3)

DATA ASSIMILATION MODELING

Proefschrift

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op maandag 6 november 2006 om 12.30 uur

door

Remus Gabriel HANEA

(4)

Samenstelling Promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. ir. A.W. Heemink Technische Universiteit Delft, promotor Dr. ir. M. Verlaan Technische Universiteit Delft, advisor Prof. dr. H. Kelder Technische Universiteit Eindhoven Prof. dr. R.M. Cooke Technische Universiteit Delft Prof. ir. N.D. van Egmond Universiteit Utrecht

Prof. dr. G.J McRae Massachusetts Institute of Technology, Cambridge, USA

ISBN

Copyright c 2006 by Remus Gabriel Hanea.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or trans-mitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the author.

Typesetting system: LA_{TEX 2}

(5)

(6)

(7)

Chapter 1 Introduction

1.1 Air pollution

Air pollution is a broad term applied to any chemical or physical (particulate matter) agent that modifies the natural characteristics of the atmosphere. The atmosphere is a complex, dynamic natural system that is essential to support life on planet earth. Stratospheric ozone depletion due to air pollution has long been recognized as a threat to the Earth’s ecosystems.

Air pollutants are classified as either directly released or formed by subsequent chemical re-actions. A direct release air pollutant is one that is emitted directly from a given source, for example the carbon dioxide or sulfur dioxide, mainly byproducts of combustion; whereas, a subsequent air pollutant is formed in the atmosphere through chemical reactions involving di-rect release pollutants.

1.2 Ozone: "

the bad, the good and the ugly

"

Since 1980, the ozone layer has got thinner, both globally and above the Netherlands. This is a result of the increased concentration of ozone-depleting substances in the atmosphere. This phenomenon is particularly evident since about 1980. Since the end of the nineties, the thick-ness of the layer has remained more or less constant at a reduced level. From 1997-2002, the ozone layer was approximately 3% thinner globally and in temperate zones such as the Nether-lands than it was prior to 1980. The reduction in the thickness of the ozone layer varies with the seasons: approximately 4% in winter and spring and approximately 2% in summer and autumn (WMO (2003)). The depletion of the ozone layer since 1980 has been particularly evident above the polar regions. In September and October, the ozone layer is on average 40-50% thinner than before 1980. Over the South Pole, depletion has stabilized since the mid-1990.

In a few recent years with cold winters in the North Pole, the ozone layer was also approxi-mately 30% thinner there. The thinning of the layer above the North Pole varies from year to year more than above the South Pole. With the depletion of the stratospheric ozone the risk to human health has increased. Therefore stratospheric ozone represents the ” good ” ozone, which keeps us protected against UV solar radiations. One of the harmful effects of UV expo-sure is the risk of skin cancer. UV expoexpo-sure has increased because of the depletion of the ozone

(12)

2 1.2 Ozone: ” the bad, the good and the ugly”

Figure 1.1: Right panel: The distribution of the ozone in the atmosphere with the altitude; left panel: the correspondence with the three types of ozone " the bad, the good and the ugly"

layer.

Ozone is also present in the troposphere. It is formed by chemical reactions involving volatile

organic compounds ( ), carbon monoxide ( ), methane ( ) and nitrogen oxides (

and , collectively referred to as ) and sunlight. Ozone is formed at different levels

in space and time. In the free troposphere, ozone is mainly formed by the reaction involving

and . Emissions come from anthropogenic sources (human activity) related to

burn-ing different kinds of fuel (combustion-fired power plants, controlled burn practices used in agriculture and forestry management and other) and natural sources (dust from natural sources, usually large areas of land with little or no vegetation, pine trees, which emit volatile organic compounds and other). The increase in emissions in recent years may be contributing to the rise in the background level of ozone in the Northern Hemisphere. Ozone in the troposphere is a greenhouse, contributing to the radiative forcing of climate and it is therefore sometimes called ” bad ” ozone [IPCC (2005)]. Under certain atmospheric conditions (sunlight, high tem-perature, little wind) photochemical smog can occur. Average hourly ozone concentrations may

rise to over 200 . Such concentrations may persist over extended periods of time over

densely populated cities, such as London, Los Angeles, Mexico City, Athens, Beijing, Hong Kong or the Ruhr Area.

(13)

Figure 1.2: The famous London smog case in 1952

ozone is not limited to the area close to where the ozone precursors are emitted. Transbound-ary fluxes transport ozone precursors over distances of hundreds of kilometers Lelieveld et al. (2002). Atmospheric chemistry transport models have been developed to understand the pro-cesses controlling the formation of ozone, to study the potential effects of ozone on ecosystems and humans, and to assess the effects of emission reductions on the ozone concentration. Because in the formation of ozone the emissions of ozone precursors plays a key role, Euro-pean countries are required to make an inventory of their national emissions of ozone precur-sors, acidifying compounds and other air pollutants every year. These emissions are reported to EMEP under the Convention on Long Range Transboundary Air Pollution (CLTRAP) and to the EU under the National Ceilings Directive [Vestreng and Goodwin (2004)]. Beside their use as input to atmospheric models for calculating concentrations and deposition of ozone, nitrogen

dioxide, s, sulfur oxides, ammonia, and other species [see e.g. Tarras´on et al. (2004)],

these national emission reports are used in assessments to ascertain if emission reductions and targets of the CLTRAP have been met.

The concentrations in the atmosphere of , , s, and other atmospheric pollutants are

(14)

4 1.3 Data assimilation

Atmospheric models and observations are usually applied separately to obtain information about ozone and air-pollutant concentrations in the boundary layer. Assimilating ozone obser-vations in the model can improve the information on the state of the atmosphere. One method to assimilate measurement data into a model simulation to get a better estimate of the real ozone concentration is data assimilation.

1.3 Data assimilation

Data assimilation (DA) is a method in which observations of the current (and possibly, past) state of a system are combined with the results from a mathematical model to produce an anal-ysis, which is considered as ’the best’ estimate of the current state of the system.

(15)

Uliasz (2003) used a Lagrangian particle dispersion model in combination with an influence function to estimate agricultural methane emissions in to local sites in New Zealand. Inverse

methods have also been applied to estimate sources and sinks of by e.g., Kaminski and

Giering (1999a,b) using an adjoint model and by Bousquet et al. (1999a,b) applying a Bayesian

three-dimensional inversion of observations. Manning et al. (2003) and O’Doherty et al.

(2004) applied a simulated annealing technique in combination with a dispersion model to

es-timate emissions of , , , and in Europe, while a Kalman filter

approach was used to estimate regional emissions of ozone depleting substances by e.g., Hart-ley and Prinn (1993) and Mahowald et al. (1997). A Kalman filter approach has also been

used by Gilliland et al. (2003) to estimate emissions of in the eastern US. There are

only a few studies in which inverse methods are applied to estimate emissions of ozone pre-cursors. Mendoza-Dominguez and Russell (2001a,b) used a four-dimensional data assimilation to estimate adjustments to emissions of ozone precursors and aerosols for the Atlanta, Georgia metropolitan area. Segers (2002) assimilated boundary layer ozone in Europe using a Kalman filter and variable emissions of ozone precursors, but did not study the emissions in detail.

1.4 Overview

The objectives and questions raised are dealt with in six chapters of this thesis. Chapter 2 in-troduces the concept of low-rank filters which one needs to use when working with large scale applications, starting with the Extended Kalman filter and finishing with the newest methods for data assimilation( the deterministic and the stochastic approach).

Chapter 3 is dealing with a very important issue for the Kalman filter applications:convergence or divergence. In this chapter we describe certain properties that a model should have, for the filter to be able to converge, and also what are the properties of the parameters that we want to be estimated

In Chapter 4 a Kalman filter approach using two low-rank algorithms (RRSQRT and ENKF) was developed to predict and estimate the ozone concentrations in combination with a large-scale atmospheric chemistry model (EUROS) above Europe. A parameter estimation for the

emissions of and VOC, the photolysis rates of and , and the deposition rate of

was performed, and the correction factors for all of the above parameters obtained. The ulti-mate goal of the data assimilation was to realize an optimal estimation of the real state of the atmosphere on the whole area and the corresponding model parameters needed for a parameter estimation tool.

Chapter 5 is a theoretical chapter where a combination of two low rank algorithms is presented and results with a real life application are described. A description of the hybrid algorithm is given and also a nonlinearity measure is introduced, allowing some insight into the performance of the different algorithms.

(16)

atmo-6 1.4 Overview

spheric chemistry transport model EUROS to estimate emissions of and s. In the

implementation of the Kalman filter with the EUROS model, model parameters (e.g. emissions) form part of the state vector. To estimate emissions a Kalman smoother technique was applied by including model parameters at previous time steps in the state vector.

(17)

Chapter 2 Overview of the low-rank filter

algorithms

2.1 Introduction

In recent years there has been a large development of data assimilation applications based on Kalman filtering in various areas of research. The classical Kalman filter algorithm is not an appropriate tool for large scale applications, due to problems encountered with the explicit treatment of the state error covariance matrix. In almost all the atmospheric chemistry trans-port models (CTMs) and in the oceanographic applications the dimensionality of the problem makes impossible the storage and the computation of such error covariance matrix. Therefore, an idea to solve the problem is to approximate the covariance matrix by a matrix of low-rank, which leads to a low-dimensional subspace of the true error space. This approach reduces the KF (Kalman filter) computational costs and in the same time solves the problem of nonlinearity imposed by the KF.

2.2 State space model

Kalman filtering represents the link between a model and measurements. This technique pro-cesses the measurements in a physically consistent way, taking into account the model. This is achieved by extending the deterministic model represented by:

(2.1) to a stochastic model:

(2.2)

where represents one time step of the numerical model and represents the model error,

which is assumed to be a stochastic perturbation with zero mean and covariance matrix .

All the available data for time are stored in a vector . The ”true” data or the ” true”

values measured without errors are related to the true state according to a linear observation model:

(18)

8 2.3 Extended Kalman filter for nonlinear systems

Through the observation model operator H, a forecast for the observed data locations can be made from the forecast of the state. Uncertainties in the measurements need to be specified as

well. Therefore, the vector from equation (2.3) is expanded as follows:

(2.4)

where is the observation error. It consists of the measurements error due to imperfect

measurements and the representation error caused by the discretization of the dynamics.

is a random vector which is assumed to be of zero mean and covariance matrix and

un-correlated with the model error .

2.3 Extended Kalman filter for nonlinear systems

For linear dynamics and for linear observation operators the classical KF approach represents the minimum variance and maximum likelihood estimator under the assumptions of model and observation errors to be Gaussian. Trying to model a real physical system, one hardly can find linear dynamics. The EKF is a first-order extention of the KF for non-linear models given by the equations (2.2) and (2.3). A detailed derivation of the equations of EKF is given in Jazwinski (1970). The only difference is that in EKF due to the nonlinearity of the dynamical operator a first-order linearization of the model operator around the estimate at the certain time is used. The algorithm can be summarized as: The first step, initialization step:

(2.5) The second step, the forecast step, defines the evolution of the distribution of the true state:

(2.6) (2.7)

(2.8) where:

(2.9) represents the tangent linear model. In general the applicability of the forecast model strongly depends on the assumption of week nonlinearity of the model.

The third step, the analysis step, in the filter algorithm is the analysis of the available data. If measurements are available at a certain time the mean and the covariance are replaced by analyzed equivalents given the new information:

(19)

where the Kalman gain is calculated by:

(2.11) These are the general equations for conventional Kalman filtering for nonlinear systems. There are some disadvantages applying this approach for real life applications.

The fact that the forecast step of EKF is based on a first order linearization, in case on nonlin-ear systems the neglection of higher order terms can lead to instabilities in the filter algorithm [Evensen (1994)]. For the linear case the KF leads to the optimal minimum variance estimate if

the covariance matrices and as well as the initial condition are completely

described. Also, the estimate is the maximum likelihood estimate. For nonlinear dynamics the EKF approach will always give an approximation of the optimal estimate. The computational costs required to store and compute the error covariance matrices in case of large scale appli-cations in atmospheric science or oceanography makes EKF not a feasible approach. A very important drawback of the EKF is represented by the fact that a tangent linear model is required. In case of large scale application this is hard to construct and translate into computer code.

2.3.1 Low-rank extended Kalman filters

These problems, encountered when one tries to implement an EKF algorithm for realistic mod-els with large state vectors lead to the development of a number of low-rank approximation algorithms. All these algorithms use a low-rank representation of the covariance matrix either by a random ensemble or by an explicit low-rank approximation of the matrix. Therefore, the filter analyses operate only in a low-dimensional subspace, denoted the error subspace, which approximates the full error space. One can pint point two different directions to achieve the approximation of the covariance matrix: factorized filters and ensemble filters.

Factorized filters

Factorized filters are based on the concept of a factorized covariance matrix. Therefore, we will consider:

(2.12) (2.13) (2.14) the factorization for the covariance of the true state, for the covariance of the forecast and the

observation error. Therefore, the covariance matrix is written as a product of a rectangular

matrix square root and its transpose:

P S S

(20)

This representation is not unique; can also be represented as:

(2.16)

where is an arbitrary orthogonal transformation. A matrix is introduced

( ) for the mapping of the forecast covariance square root to the observation space. In this

Chapter we will make the following notation . With these notation the forecast of

mean and covariance can be computed as:

(2.17) (2.18)

or As one can see the introduction of the forecast error leads to an extension

of the square root with the columns of . But, at the same time every column introduces a new

direction for the uncertainty in the state vector. To stop the number of modes to grow with every time step one needs to include truncation mechanisms. The equations for the analysis of the filter can be computed with the minimal variance gain. The analysis equations for a minimal variance gain are:

(2.19) (2.20)

(2.21) Or, we can rewrite:

(2.22) Although the degree of freedom in the state vector is very large, the errors in the state are described very well by a limited number of directions. It does not matter if these directions are represented by singular vectors, modes , empirical orthogonal functions (EOFs), the underlying equations and implementation remain the same.

Singular Evolutive Extended Kalman filter (SEEK) and Singular Evolutive Interpolated Kalman filter (SEIK)

(21)

by a limited number of EOFs. These EOFs are the eigenvectors of a sample covariance matrix , computed over a large number of state vectors obtained with the deterministic

model. If one denotes obtains a covariance square root formulation. For SEEK

filter the noise is projected onto . In case of SEEK filter the mean and covariance base are

propagated separately by a tangent linear model and in case of SEIK, for the propagation of the

covariance matrix an ensemble of states from the columns of are formed.

Error subspace Statistical Estimation (ESSE)

The Error subspace Statistical Estimation has been introduced by Lermusiaux and Robinson (1999) and represents a framework for data assimilation techniques based on low-rank approx-imation of the true error covariance matrix. Taking into account the problems which appear in case of ocean-atmosphere applications lead to a definition of the filter problems in terms of flexible sized error subspaces that spans and follows the scales and physical processes where the most important errors occur.

Depending on the application (filtering, forecast, smoothing parameter estimation) one can choose the proper technique for low-rank filtering. The main idea is to approximate the

er-ror covariance matrix by a factorization where is a base for the error subspace

and contains the eigenvalues according to some norm. One can reduces the rank of the error

subspace by using singular value decomposition. The forecast step is based on the propagation of an ensemble of states chosen random or preconditioned by the mean and the error subspace.

2.3.2 Ensemble Kalman filters

Ensemble data assimilation methods assimilate observations using state-space estimation meth-ods and low-rank representations of forecast and analysis error covariances. The main idea is to transform the forecast ensemble into an analysis ensemble with the correct statistics. One can obtain this stochastically by treating observations as random variables, or deterministically, by requiring that the update analysis perturbations satisfy the Kalman filter analysis error covari-ance equations.

Stochastic approach

Where the factorized type of filters are based on factorization of the covariance matrix, the En-semble Kalman filter is based on the convergence of large numbers. The EnEn-semble Kalman filter was introduced by Evensen (1994) and has been successfully used in many applications [Evensen and van Leeuwen (1996); Houtekamer and Mitchell (1998)]. This Monte Carlo ap-proach is based on the representation of the probability density of the state estimate in an

en-semble of possible states, . Each ensemble member is assumed to be a single

(22)

generated system states (also called modes).

In the first step of this algorithm an ensemble of N states is generated to represent the

uncertainty in . In the second step, the forecast step, the stochastic model propagates the

distribution of the true state:

(2.23) When the measurements become available the mean and the covariance are replaced with equiv-alent ones in the analysis step using:

(2.24)

where the ensemble of state vectors is generated with the realizations and of the

noise processes and , respectively. The advantages of this algorithm are that

is positive definite and that the linear tangent model is not required any more because the en-sembles are propagated through the model using the original operator. Another advantage is the fact that there is no need for rank reduction, since the ensemble size does not grow due to the introduction of system noise. And also, in the final implementation of the algorithm, need not be computed [Evensen (1994)]. For most practical problems the forecast equation is computationally dominant [Dee (1991)]. As a result, the computational effort required for the ENKF is approximately N model simulations. The standard deviation of the errors in the state estimate is of a statistical nature and converge very slowly with the number of ensembles. This is a disadvantage of a Monte Carlo approach.

In the ENKF type of algorithms , if the analysis error covariances systematically underestimate ensemble mean analysis error, observations will be underweighted, and filter divergence may occur. In order to avoid this problem methods for increasing the ensemble variance are neces-sary.

Many approaches were developed in order to solve this tendency toward filter divergence. One approach is to localize background error covariances by applying a Schur product with a cor-relation function, as discussed in Houtekamer and Mitchell (2001). The product of these two covariance models reduces spurious noise in the covariances and the resulting tendency toward introducing unrealistically large analysis increments for from the observation.

Another stochastical approach to solve the problem of filter divergence is to use the”double”

ENKF (DENKF) [Houtekamer and Mitchell (1998)], where ensemble members are split in two

(23)

This can help prevent the feedback cycle which leads gradually to very small forecast error co-variances.

A hybrid filter which calculates the covariances as combinations of covariances from the en-semble and from a stationary model like 3DVAR was introduced by Hamill and Snyder (2000). Deterministic approach

There are two fundamental problems associated with the use of the ENKF as written above. First, ensemble size is limited by the computational cost of applying the forecast model to each ensemble member. In the same time, small ensembles have few degrees of freedom available to represent errors and suffer from sampling errors that will further degrades forecast error co-variance representation. Sampling errors leads to loss of accuracy and underestimation of error covariances. This problem can progressively worsen, resulting in useless ensemble or forecasts, which can cause filter divergence.

If one already has a problem that suffers from sampling errors, the addition of perturbed obser-vations adds an additional source of sampling errors related to the estimation of the observation error covariances. In addition to reducing the accuracy of the analysis error covariance estimate, this extra source of sampling error increase the probability that the analysis error covariance will be underestimated. Therefore, an ensemble data assimilation scheme which uses no perturbed observations are expected to perform better than the ones with perturbed observations for a cer-tain type of applications.

Using the ”Potter method” for the square root Kalman filters [Bierman (1977)] one can obtain equation (2.21) and rewrite as:

(2.25) (2.26) (2.27)

where we defined the matrix .

Then the analysis perturbation ensemble is calculated from:

(2.28)

where: and is an arbitrary orthogonal matrix.

Therefore, one can say that the update ensemble is a linear combination of the columns of

and is obtained by inverting the matrix and computing a matrix square root

of the matrix When observation errors are uncorrelated, observations

can be assimilated one at the time or serially [Houtekamer and Mitchell (2001); Bishop et al.

(2001)]. For a single observation, , is a column vector, and the innovation is a scalar.

In this case, a matrix square root of can be computed as:

(24)

and solving for the scalar , which gives:

(2.30)

The analysis ensemble update for is:

(2.31) At observation locations, the analysis error ensemble is related to the forecast error ensemble by:

(2.32)

The scalar factor has an absolute value less than or equal to one and it is

positive when the plus sign is chosen in the definition of .

In Whitaker and Hamill (2002) the analysis perturbation ensemble is found from:

(2.33)

where the matrix is a solution of the nonlinear equation:

(2.34) In case of a single observations a solution of the above equation is:

(2.35)

where the plus sign is chosen in the definition of . The corresponding analysis perturbation

ensemble update:

(2.36)

which is identical with equation (2.31). This algorithm is calledEnsemble Square Root Kalman

filter (ENSRF).

Another method of computing the update analysis ensemble is to use the Sherman-Morrison-Woodbury identity to show that:

(2.37)

The matrix on the right-hand inside of the above equation is practical to compute when

the inverse observation error covariance matrix is available. This approach avoid inverting

the matrix and it is used in the Ensemble Transform Kalman filter (ETKF)

(25)

Anderson (2001) proposed a method denotedEnsemble Adjustment Kalman filter (EAKF)

where the analysis is computed without adding perturbations to the observations. If observations are not perturbed in the ENKF this still gives the correct mean of the analyzed ensemble, but results into a lower variance as explained by Burgers and Evensen (1998). This in EAKF is accounted for by deriving a linear operator which replaces the traditional gain matrix and results in an update ensemble which is consistent with theory. The algorithm is based on the increasing background error covariances by inflating the deviation of forecasted members with respect to their mean by a small amount. Before the first observation is assimilated in a new cycle, forecast ensemble deviation from the mean are inflated by an amount , slightly greater

than 1.0 ( ).

(2.38) The purpose of this replacement is to account for a slight under representation of variance due to the use of a small ensemble. The fact that the square root representation of the error covariance

matrix is not unique, beginning with the same forecast error covariance, the ETKF, EAKF

and ENSRF methods produce different analysis ensembles that span the same state-subspace and have the same covariances.

Reduced Rank Square Root Kalman Filter

The RRSQRT filter was developed by Verlaan and Heemink (1995) and was successfully used in different types of applications: water level measurements in a shallow water model and to a number of hydro-dynamical models [Canizares (1999)].

In the formulation of this algorithm, the covariance matrix is expressed in a limited number of orthogonal modes, which are re-orthogonalized and truncated to a fixed number (initially choose) during each time step. An important step in the RRSQRT algorithm is the reduction of the covariance square root. After the introduction of the noise system the number of modes has

grown from to , where is the number of columns in . The reduction step reduces

the number back to using a singular vector decomposition approach:

(2.39) (2.40)

where matrix contains the eigenvectors of corresponding with the largest

eigenval-ues and the diagonal matrix contains the eigenvalues in descending order. The new formed

matrix is an approximation of maintaining the largest singular vectors.

(26)

(27)

Chapter 3 Convergence and divergence in the

Kalman filter approach

1

3.1 Introduction

In the past years applications of Kalman filter theory were noted and successfull in many areas of research, starting with the meteorological applications [Ghil and Malanotte-Rizzoli (1991);Courtier et al. (1993)], tidal models [Heemink and Kloosterhuis (1990)], nonlinear shallow-water storm-surge models [Verlaan and Heemink (1996)] and in atmospheric chem-istry and transport modeling has received attention as well [e.g., Van Loon and A.W.Heemink (1997); Segers et al. (2000); Elbern et al. (1997); Elbern et al. (2000);Wang et al. (2001)]. As these problems start to become more complicated form the physical point of view, or to become a large scale application some of the classical Kalman filter approaches started to suffer of di-vergence.

An important task is to point the causes for this divergence problem and try to eliminated as accurate as possible.

As the Kalman filter theory was originally written (Kalman) an important property of the model is to be linear. But, in real life applications rarely one can use a linear model to describe very complicated physics. Or, if one considers the recursive parameter estimation problem it ends up with a nonlinear problem. The next step in solving this nonlinearity of the model was the Extended Kalman filter which is an approximate filter based on first-order linearization. It was used for joint parameter and state estimation problems. The essential difficulty in case of approximation techniques, as EKF, is to establish convergence and Ljung (1979) gives an overview of the causes of divergence. The main cause is that the EKF is not able to handle the nonlinearities very well. Also Verhaegen and Van Dooren (1986) showed that the numerical

er-rors in EKF leads to divergence. Evensen and van Leeuwen (1996) introduces a new approach

for the Kalman filter theory, the Ensemble Kalman filter. This technique solves successfully the problem of nonlinearity of the model,without the first-order linearization, which was the cause of the numerical errors. Still, often successful, the ENKF gives an approximation, even

1_{This Chapter it is based on the article ”Convergence and divergence in the Kalman filter approach” submitted}

to MWR in November 2006

(28)

18 3.2 The theory behind

for , where is the number of members in the ensemble.

Another reason for filter divergence is represented by the modeling inaccuracies as Lewis (1986) shows.. Modeling inaccuracies, including noise statistics, can lead to nonconvergence and un-bounded actual estimation errors and can lead to biased estimates.

But, if one can model perfect the physics of a large scale application, there will be problems with working with huge and ill conditioned error covariance matrices. One of the solution of this problem consists in low rank Kalman filter algorithms, which are based on some ap-proximations of the error covariances matrices, by orthogonally projecting them on a smaller space. The goal of this paper is to show that this kind of approximations of the error covariance matrices, includes errors in the sequential algorithm that can grow and leads in some cases to divergence. The key parameters of this divergence problem and the ways that one can tune it are shown here. In Section 3.2 the basic presumptions made in the system and control theory for state and parameter estimation are presented. In Section 3.3 some simple examples are shown to picture the possible inaccuracies that can be encountered. An analysis of the possible errors in the approximation of the error covariance matrices are presented in Section 3.4. In Section 3.5 the same simple examples presented in Section 3.4 are generalized for a more complicated twin experiment based on a 2D advection diffusion equation. In Section 3.5 the results of this twin experiment are presented and the conclusions.

3.2 The theory behind

Kalman filter is an ideal tool when all the dynamics and noise statistics are exactly known. If this is not the case some serious divergence problem can occur, because of the differences between the theoretical behavior of the filter and its actual behavior. We define two concepts. The ”truth” is the actual physical system that we are interested in. The model is the mathematical description, which is used for the filter design stage.

Let the model be given by:

(3.1) (3.2)

with and white and uncorrelated with each other and with . Let ,

and

The Kalman filter designed on the basis of this model is given by:

Time update:

(29)

Measurement update:

(3.5) (3.6) (3.7) The a priori error covariance is defined as:

(3.8) Now suppose that in reality the ”truth” is described not by equations (3.1) and (3.2) but by the

state with the following dynamics:

(3.9) (3.10)

with and white and uncorrelated with each other and with . Let ,

and

One need to develop a theoretical background for analyzing mistmatch between the model and the ”truth”. We define the a priori and a posteriori between the estimate generated by (3.3)-(3.7) and the state of the actual ”truth” as:

(3.11) (3.12)

If the filter is designed and implemented correctly the actual errors and should be

accept-ably small. Therefore, is interesting to calculate the expressions for the actual error covariances

and .

First consider the time update:

(3.13) We denote the model mismatch system matrix:

(3.14) and we can rewrite equation (3.13) as:

(30)

20 3.2 The theory behind

We define the augmented state vector as:

(3.16)

and we obtain the following equality:

(3.17)

We define the correlation matrix:

(3.18)

From equation (3.17) we obtain:

(3.19)

where we have set and for simplicity. Multiplying this equation the result is the

time update equations:

(3.20) (3.21)

(3.22)

Note that is the ”truth” state correlation matrix .

If , and , so there is no model mistmatch, then equation (3.22) is

identical with equation (3.3). The extra terms in Eq. 3.22 represent the increase in the error

covariance due to the modeling inaccuracies. The only quantity of interest is , which is the

correlation matrix of the actual error . As one can see, in the presence of modeling errors the

actual error correlation depends on the ”truth” state correlation. Now we will derive the measurement update equations:

(3.23)

One can observe that in the above equation was used and not the model output . This is

(31)

therefore the actual ”truth” output. We can rewrite:

(3.24) where we denoted:

(3.25) By defining an augmented state vector we can rewrite equation(3.24) as:

(3.26) If we take the correlation in both sides we obtain:

(3.27) Therefore, the measurement update equations are:

(3.28) (3.29)

(3.30)

Note that if and so that are no modeling errors, then the update equation (3.30)

for the actual error correlation is equal to update equation (3.6) for the predicted covariance. We

are not interested in and , but we need these quantities if we want to find the quantity of

.

Since initially the uncertainty in the estimate is the same as the uncertainty in the ”truth” state, we would initialize the recursion by:

(3.31) By combining equation(3.15) and equation (3.24) we obtain the error system for the actual error:

(32)

22 3.3 Illustrative examples

Hence:

(3.33) There are several comments to be made about this equation:

It is driven by the ”truth” noise process.

If the model mismatch matrices and are not both zero, the error system is driven

by the ”truth” state . Unless is bounded, the error will grow. Therefore, in

the presence of modeling inaccuracies, the ”truth” must in most cases be stable to avoid causing an unbounded estimation error.

The effect of poorly known deterministic input is also clearly seen in equation (3.33),

for the error system is driven by the term . If this difference is not zero,

unbounded errors can again occur in .

One can conclude that even if is stable for large , which is guaranteed by the

steady-state properties of the Kalman filter, modeling inaccuracies can lead to nonconvergence and unbounded actual estimation errors. Modeling inaccuracies can lead to biased estimates.

3.3 Illustrative examples

We will show three hypothesis in which one can get into divergence problems and we will present the conditions (sufficient and necessary) for the system and for the model characteristics in case of a divergent or convergent filter.

Example1

To estimate the value of a constant scalar unknown from the measurements corrupted by

additive white noise, the system and the measurement model are:

(3.34) (3.35)

where and and is white. Let

Before any measurements are taken one can determine the error covariance off-line.

(33)

– Time update (3.36) – Measurement update (3.37) – Time update (3.38) – Measurement update (3.39) .. . – Time update (3.40) – Measurement update (3.41)

To find the estimate incorporating the a-priori information on and all measurements

(34)

24 3.3 Illustrative examples 0 50 100 150 0 0.5 1 k pk

Figure 3.1: Behavior of the predicted error covariances

From equation (3.40):

(3.46) Therefore, using equations (3.43) and (3.46), equation (3.45) becomes:

(3.47)

The predicted error covariance is plotted in Figure 3.1 for and ; it shows the

anticipated decrease to zero of the error as more and more measurements are made. The Kalman gain also goes to zero with ; this means that as more measurements are made, the estimate theoretically becomes reliable so that each subsequent measurement has less and less effect. Now, let the actual ”truth”, which is supposedly modeled by (3.1) and (3.2), be given by:

(35)

where , , . Thus, the actual unknown is a Wiener process, which we have incorrectly modeled as a random constant. To find the resulting actual error correlation for:

, we need to use the propagation formula:

(3.50)

(3.51)

Since and , we do not need to compute , but simply use the

recursion for the actual error correlation .

(36)

26 3.3 Illustrative examples

– Time update

(3.54) – Measurement update

(3.55) In general we can write:

(3.56)

If we take , , the actual error correlation is plotted in the same picture

with the predicted covariance . As one can see, the behavior of the actual error correlation is

quite different from the behavior of the predicted covariance. To find the limit of the behavior

of one can use:

(3.57) (3.58) and obtain: (3.59) (3.60) (3.61) (3.62) (3.63)

One can observe that while for large everything seems well because of the eq. 3.41 with

(37)

0 50 100 150 0 0.5 1 k p a and sk a k sa_k pa k K_k

Figure 3.2: Behavior of the predicted error covariances

3.4 Impact of suboptimal filter algorithms

In the first example the errors were due to the wrong modeling of the actual unknown Wiener process, which we have incorrectly modeled as a random constant. This example will represent the situation when we model correctly an actual constant unknown, and the difficulties appear from the approximation of the covariances. In this example, let the actual ”truth” be given by:

(3.64) (3.65)

where , . Thus, the actual unknown is a random constant. And the

model is given by:

(3.66) (3.67)

where , . Thus, we are in the situation where the model is perfect

(38)

28 3.4 Impact of suboptimal filter algorithms

consider the computed Kalman gain and computed error covariance matrix:

(3.68) (3.69) (3.70) Now, if we return to the ”truth” and write down the propagation equations for the real error, using the computed Kalman gain:

(3.71) (3.72) Equation 3.72 can be obtain from 3.30 with

3.4.1 Underestimation of the error covariance matrix

In case of large scale applications the full Kalman filter is not practical to be implemented, and the use of suboptimal filters is a way to bypass the problems. The suboptimal filters are based on the idea of approximating the full covariance matrix, projecting the propagation equations onto a smaller supspace. In case of the RRSQRT filter (discussed in Chapter 2) one can prove

that filter error covariance P SS systematically underestimates the covariance. Moreover

the following lemma can be proven [Verlaan (1998)]. Lemma

Consider:

P = covariance of the estimation error of the full model

P = reduced rank approximation of P

P = actual covariance if a suboptimal reduced rank filter is used in the full model Then we have the following result:

(3.73)

Using the small example that we have presented in the previous section we systematically

un-derestimating the covariance (we chose to be negative always). The results can be seen in

(39)

0 10 20 30 40 50 60 70 80 90 100 ï1 0 1 2 3 4 5 6 k (time step) Pk , S k , O k calculated erros true errors optimal errors 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 k (time step) Pk , S k calculated erros true errors optimal errors

Figure 3.3: Behavior of the predicted error covariances for the underestimating (deterministic case against stochastic case)

right panel we have the stochastic case, which allows model error. In the case of weak constrain one can observe the same result as in the Lemma and also the fact that the filter is stable.

The other case is when the approximation randomly underestimates or overestimates the error covariance matrix. The results obtain are plotted in Figure 3.4 with the left panel we have the deterministic case without model errors, and on the right panel we have the stochastic case, which allows model error. In both cases presented above when trying to estimate an constant parameter using the deterministic model and approximating the error covariance matrix lead to divergence. Introducing stochastic forcing helps in both cases.

3.5 Illustrative examples (twin experiment)

Before using the algorithm for real-life atmospheric chemistry data assimilations problems, first we define a twin experiment .

3.5.1 The model

Consider the 2D advection diffusion equation:

with a square domain and zero initial conditions. Here, is the concentration,

(40)

30 3.5 Illustrative examples (twin experiment) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 k (time step) Pk , S k , O k calculated erros true errors optimal errors 0 10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 k (time step) Pk , S k , O k calculated erros true errors optimal errors

Figure 3.4: Behavior of the predicted error covariances for the random under and over estimation (deterministic case against stochastic case)

backward Lagrangian scheme is used to discretized the equation on a grid.

In addition, we consider a simple chemistry which describes the conversion of to

using the following chemical reactions:

where the rates and are considerated to be constants. The

second reaction is used to represent cloud chemistry and oxidation path-ways. We have also

considered a new coefficient which represents a constant deposition rate for

component. The change in concentration fields can be describe as follows:

Having two processes: transport and chemistry, both in the same time step, we use splitting operator technique to separate the processes and solve them separately.

In order to compare the algorithms to other some twin experiments were carried out. A reference

solution was generated by inserting constant emissions at five grid cells The increase of

(41)

The final product is . Therefore the observations were generated from simulated truth

concentrations of , which were computed by adding fluctuations to the mean emissions

according to:

where is an independent Gaussian white noise processes with and var

.

The index refers to source locations

. Negative emissions are truncated and white observational noise

(with variance ) is added to the truth concentrations.

The state vector here consists of 1805 compounds :

concen-trations of , concentrations of and 5 points of pollutant emission.

3.5.2 Uncertainties in emissions

In our model the emissions of have been modelled uncertainly. With the deterministic

value of the emission in a grid cell, the stochastic emissions are modelled according to:

We treat the emissions as a model input and we define our model in the following:

(3.74)

The noise on the model and on the measurements are fixed for the simulated data. We assume we do not know them and we use the filtering and smoothing procedures with several different noises in order to study the smoother behavior.

We know that is already well estimated using ENKF, since it is the pollutant that we

mea-sure. Therefore, we predict that the best improvement brought by the smoother algorithms will

be done on the estimates of and on the emissions. All the experiments are made in the

same initial conditions, using ensemble members and a set of measurement points.

Case

The reference solution is considered to be perfect. The parameter is considered perfect:

(42)

32 3.5 Illustrative examples (twin experiment) 0 50 100 150 3 4 5 6 7 8 9 10 11 12 so2 concentration truth std std filter 0 50 100 150 ï10 ï5 0 5 10 15 20 25 second emiss truth std std filter r=0.01 50 ens 0 50 100 150 ï2.5 ï2 ï1.5 ï1 ï0.5 0 0.5 1 1.5 2 2.5

The inovations of the filter

t

Innov

Innovations Minus theoretical standard error Theoretical standard error

Figure 3.5: Illustration of deterministic case; top left panel: concentrations in one grid cell; top right panel: second emission point and second row center: the innovations for the corresponding emission

Figure 3.5 shows that the filter cannot estimate the reference constant solution properly. In addition, the standard deviation is very small, the filter solution is convergent, but to the wrong

estimate. In Figure 3.5 we are presenting the estimate for the concentration in one grid

cell in the same conditions as in Case . Another emission is plotted in the top right panel to show that the behavior saw in the previous plot stands also for another locations for emissions. Adding uncertainty in the observations makes the filter to be more stable and one can observe that the innovations are in between the theoretical statistic. Studying the behavior of the filter in order to estimate the constant emission parameter using the deterministic model, we may conclude the following:

Starting with different initial conditions the filter estimate converges (as it is shown by standard deviation plot) but to the wrong solution.

(43)

0 10 20 30 40 50 60 70 80 1 2 3 4 5 6 7 so2 concentration truth std std filter 0 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 second emiss truth std std filter 50 ens q=0.01 0 10 20 30 40 50 60 70 80 ï2 ï1.5 ï1 ï0.5 0 0.5

t

Innov

Figure 3.6: The illustration of the stochastic case. Top left panel: concentrations in one grid cell; top right panel: second emission point and second row center: the innovations for the

corresponding emission. Using EnKF using ensembles, ,

Increasing the noise in the observations, the deterministic filter provide better results. In order to avoid the divergence of the filter, we consider the model to be stochastic.

Case

The reference solution is generated in the same way:

The model is stochastic adding the white noise in the input emissions (q=0.01):

(44)

34 3.6 Conclusions 0 50 100 150 200 250 ï2 0 2 4 6 8 10 12 14 16 18 second emiss truth std std filter 0 50 100 150 200 250 ï2.5 ï2 ï1.5 ï1 ï0.5 0 0.5 1 1.5 2 2.5

t

Innov

Figure 3.7: Time evolution of emissions points using EnKF using ensembles, , , stochastic model.

(Figure 3.7). In comparison if we look at the second emission estimate in Figure 3.6 and Figure 3.7 increasing the noise leads to a faster convergence rate to the real value. The filter recovers faster. For the last case we consider the stochastic reality and model.

Case

The reference solution is considered to be stochastic:

where variances for the white noise processes (observations and emissions) are . The model is considered to be stochastic as well:

For the last case only the third emission is presented in Figure 3.8. It can be observed that the fact that both the model and the truth is stochastic does not bring too much improvement in filter performance.

3.6 Conclusions

(45)

0 5 10 15 20 25 30 35 40 45 50 20 25 30 35 40 45

50 third emiss/ alpha=0.9

truth std std filter

Figure 3.8: True and computed concentration fields using EnKF after time steps using ensembles.

(46)

(47)

Chapter 4 Kalman filter with a large scale

application

1

A Kalman filter coupled to the atmospheric chemistry transport model, EUROS, has been used to estimate the ozone concentrations in the boundary layer above Europe. Two Kalman filter algorithms, the Reduced Rank Square Root (RRSQRT) and the Ensemble Kalman filter (ENKF), were implemented in this study. Both required, in general, a large number of EU-ROS model simulations for an assimilation. The observations consisted of hourly ozone data in a set of 135 ground-based stations in Europe for the period, June 1996. Half of these sta-tions were used for the assimilation and the other half only for validation of the results. The combination between data assimilation (Kalman filter) and the atmospheric chemistry transport model, EUROS, gave more accurate results for boundary layer ozone than the EUROS model or measurements used separately. The average difference between assimilated and measured

ozone concentrations decreased from 27.4 to 20.5 for the average of the stations used

for validation in Europe. Both algorithms tend to converge to about the same accuracy, with an increasing number of EUROS model runs. About 10-20 EUROS model calculations were found sufficient for a good assimilation. The results are supported by a number of simulations that also reveal a local character for the assimilation process.

4.1 Introduction

Elevated concentrations of ozone in the boundary layer can cause adverse affects to humans and ecosystems EC (1996). Ozone in the boundary layer is formed by chemical reactions of

the ozone precursors, nitrogen oxides ( ), volatile organic compounds (VOCs), and carbon

monoxide and methane under the influence of sunlight. The impact of ozone is not limited to the area close to where the ozone precursors are emitted. Transboundary fluxes transport these precursors over distances of hundreds of kilometers [Lelieveld et al. (2002)]. Atmospheric chemistry transport models have been developed to understand the processes controlling the formation of ozone, to study the potential effects of ozone on ecosystems and humans, and to

1_{This chapter is based on the article published in JGR: Hanea et al. 2004 ” Data assimilation of ground-level}

ozone in Europe with a Kalman filter and chemistry transport model, Journal of Geophysical Research, Vol 109, D10302, doi:10.1029/2003JDD004283

(48)

38 4.1 Introduction

assess the effects of emission reductions on the ozone concentration. Atmospheric models and observations are usually applied separately to obtain information about ozone and air-pollutant concentrations in the boundary layer. Assimilating ozone observations in the model can im-prove the information on the state of the atmosphere. One method to assimilate measurement data into a model simulation to get a better estimate of the real ozone concentration is data assimilation.

The most well-known application of this method has been in weather forecasting problems, in which data assimilation was applied in real life applications for the first time [Ghil and Malanotte-Rizzoli (1991); Courtier et al. (1993)]. The method has already proved to be useful in other fields of application like tidal models [Heemink and Kloosterhuis (1990)], or nonlinear shallow-water storm-surge models [Verlaan and Heemink (1996)]. In the last few years, data assimilation in atmospheric chemistry and transport modeling has received attention as well [e.g., Van Loon and A.W.Heemink (1997); Segers et al. (2000); Elbern et al. (1997); Elbern et al. (2000); Wang et al. (2001)]. The variational and sequential methods rank as two of the most well-known methods for applying data assimilation. Variational methods [Talagrand and Courtier (1987)] are based on the minimization of the cost function for the difference between the model and the measurements. The problem of minimizing this function is efficiently imple-mented by solving the adjoint problem. In Elbern et al. (1997) variational data assimilation is implemented using a 4D-var smoother technique for a tropospheric chemical transport model. The disadvantage of this method is the necessity of writing an adjoint model. For large-scale atmospheric transport problems implementation of an adjoint model is a very complicated task. From the beginning of 1980s, the same time as the variational methods were introduced, se-quential methods started to attract more and more attention [Ghil et al. (1981)]. At first, Kalman filter techniques were used especially in oceanographic applications [Heemink and Kloosterhuis (1990); Evensen (1994)]. The successful implementation of this method led to the introduction of similar techniques in related fields such as coastal hydrodynamics [Verlaan (1998)] and air pollution [e.g., Zhang and Heemink (1995); Segers (2002)].

Nevertheless, in the case of atmospheric transport models with a very large number of elements in the state vector, a filter approach with a full covariance matrix will be impossible to use. Sev-eral specific formulations of the Kalman filter technique were constructed to get pass this prob-lem. These are called low-rank Kalman filter algorithms, and in this research are represented by the RRSQRT and ENKF. The performance of these two algorithms was tested in Canizares (1999) and compared in a twin experiment using a two-dimensional shallow-water equation model. In the move from the twin-experiment stage in Zhang et al. (1999), the RRSQRT-filter was successfully used in real-life application to assimilate data into an atmospheric transport

model in order to evaluate the budget in Europe. The model implemented for this

re-search is a linear one and has no chemistry involved. In Segers (2002) the RRSQRT and ENKF algorithms were used to estimate and predict ozone concentrations in an area of northwest Eu-rope (England and Wales). The algorithm performance was tested on a small-scale atmospheric

(49)

six days and the measurements taken from 11 rural sites. Important model errors in the case of ozone prediction are due to uncertainties in the emissions of pollutants. To take this into

account, the , VOC and CO emissions in Segers (2002) were considered as elements of the

state vector and estimated by the Kalman filter technique. The results were promising and the data assimilation methods proved to be able to estimate the ozone concentrations in the case of a real-life small-scale atmospheric model. However, the methodology in this research field still has to be applied and modified for a real-life large-scale atmospheric model.

In our study a Kalman filter approach using the same two low-rank algorithms (RRSQRT and ENKF) was developed to predict and estimate the ozone concentrations in combination with a large-scale atmospheric chemistry model (EUROS) above Europe. The EUROS model grid

consists of cells. The assimilation period was June 1996 and 135 ground-based

ob-servational sites spread all over Europe were used. A parameter estimation for the emissions

of and VOC, the photolysis rates of and , and the deposition rate of was

per-formed, and the correction factors for all of the above parameters obtained. The emissions parameters and the deposition rates were also spatially defined for four large regions in Europe, chosen on the basis of their importance. The influence of these parameters on the performance of the assimilation was studied. The ultimate goal of the data assimilation was to realize an optimal estimation of the real state of the atmosphere on the whole area and the corresponding model parameters needed for a parameter estimation tool.

The Kalman filtering technique is discussed in section 4.2, as well as the implementation of the two algorithms. A short description of the EUROS model is given in section 4.3, followed by the definition of the statistical environment, and the observations used in the simulations. The results of the data assimilation calculations are presented and discussed in section 4.4 and the conclusions in section 4.5.

4.2 Kalman filtering

4.2.1 State space model

Kalman filtering represents the link between a model and measurements. This technique pro-cesses the measurements in a physically consistent way, taking into account the model. This is achieved by extending the deterministic model represented by:

(4.1) to a stochastic model

(4.2)

where represents one time step of the numerical model and is the white-noise

process: . The covariance matrix is supposed to be of low rank so

that it can be factorized as where is a matrix with a small

(50)

40 4.2 Kalman filtering

specified.

The state space representation in this research was based on an augmented state vector described

in detail in Section 4.4. All the available data for time are stored in a vector . The ”true”

data or the ” true” values measured without errors are related to the true state according to a linear observation model:

(4.3) Through the observation model operator H, a forecast for the observed data locations can be made from the forecast of the state. Uncertainties in the measurements need to be specified as

well. Therefore, the vector from equation (4.3) is expanded as follows:

(4.4)

where the observations operator assigns the ozone concentration for the observation

height of a grid cell to an observation and is the observational noise process:

.

4.2.2 Extended Kalman filter for nonlinear systems

The idea is to combine the measurements modeled by equation (4.4) with the information pro-vided by the atmospheric model in order to obtain an optimal estimation of the state of the system. The first step in a Kalman filter algorithm is to specify the initial distribution of the true state

(4.5) The second step, the forecast step, defines the evolution of the distribution of the true state:

(4.6)

(4.7) where:

(51)

analyzed equivalents given the new information:

(4.9) (4.10)

where the Kalman gain is calculated by:

(4.11)

These are the general equations for conventional Kalman filtering for nonlinear systems. Using these equations in an atmospheric chemistry transport model with nonlinear processes and with

a large grid involves certain problems. is very big (high computation time) and is ill

conditioned [Heemink et al. (2001)], that is,there is a large range of eigenvalues. Furthermore the tangent linear model given by equation (4.8) must be used. To solve these problems one needs to use some low-rank Kalman filter has to be used suitable for data assimilation in models with large state vectors. Two algorithms are used for this in combination with the EUROS model: the Reduced Rank Square Root Kalman filter (RRSQRT) and Ensemble Kalman filter (ENKF).

4.2.3 Ensemble Kalman filter

The Ensemble Kalman filter was introduced by Evensen (1994) and has been successfully used in many applications [Evensen and van Leeuwen (1996); Houtekamer and Mitchell (1998)]. This Monte Carlo approach is based on the representation of the probability density of the state

estimate in an ensemble of possible states, . Each ensemble member is assumed

to be a single sample out of a distribution of the true state. Whenever necessary, statistical moments are approximated with sample statistics. Thus, the ensemble Kalman filter algorithm is based on a representation of the probability density of the state estimate by a finite number of N randomly generated system states (also called modes).

In the first step of this algorithm an ensemble of N states is generated to represent the

uncertainty in . In the second step, the forecast step, the stochastic model propagates the

distribution of the true state:

(52)

42 4.2 Kalman filtering

When the measurements become available the mean and the covariance are replaced with equiv-alent ones in the analysis step using:

(4.15) (4.16) (4.17)

where the ensemble of state vectors is generated with the realizations and of the

noise processes and , respectively. The advantages of this algorithm are that

is positive definite and that the linear tangent model is not required any more because the en-sembles are propagated through the model using the original operator as in equation (4.12).

And also, in the final implementation of the algorithm, need not be computed [Evensen

(1994)]. For most practical problems the forecast equation (4.12) is computationally dominant [Dee (1991)]. As a result, the computational effort required for the ENKF is approximately N model simulations. The standard deviation of the errors in the state estimate is of a statistical nature and converge very slowly with the number of ensembles. This is a disadvantage of a Monte Carlo approach. Here it should be noted that for many atmospheric data assimilation problems the analysis step is also a very time consuming part of the algorithm [Houtekamer and Mitchell (1998)].

4.2.4 Reduced Rank Square Root Kalman filter

Another approach for solving large-scale Kalman filtering problems is to approximate the full covariance matrix of the state estimate by a matrix with a reduced rank. This approach was introduced by Cohn and Todling (1995) and Verlaan and Heemink (1995), Verlaan and Heemink (1996). Algorithms based on similar ideas have been applied by Pham et al. (1998a). The

reduced-rank approach can also be formulated as an ENKF where the ensemble members

are not chosen randomly, but in the direction of the leading eigenvectors of the covariance

matrix [Verlaan and Heemink (1997)]. This algorithm can be applied and implemented in case

of the matrix is a low rank matrix. Because of the fact that few model parameters were

assumed to be uncertain and because of the dispersion in the model we can assume that is

low rank in our case. The algorithm approximates the error covariance for forecast and

analysis by:

(4.18)

where contains the eigenvectors of the covariance matrix and is the diagonal matrix

with eigenvalues on the diagonal. The square root of is denoted by:

(4.19)

So, will be a matrix because it will take the first -eigenvalues of the matrix. Now

(53)

propagation of the estimate and of the error covariance take place according to:

(4.20)

The addition of columns every time step, represented by for the noise, would increase the

computation time.

The second step in the algorithm is the reduction step: the idea is to use only the first q leading

eigenvalues and eigenvectors of the approximated covariance matrix . It is known that

has the same eigenvalues as the matrix and that is why a

singular-value decomposition is carried out:

(4.21) This method places the eigenvalues in decreasing order and can also prove that the eigenvectors

of are given by the columns of the matrix . The actual computation

is then given by the eigenvalue decomposition of equation (4.21) followed by a reduction like: (4.22)

Here, denotes the truncation of the number of columns to .

The third step is the measurement update when observations are available. Several of the mea-surement update equations known in literature do not depend on the specific type of square root. Here is presented the scalar update of Potter [Maybeck (1979)]:

K (4.23)

K

Independent measurements can be processed one at a time. The most time-consuming part of the RRSQRT algorithm is the propagation of the modes (equation (4.20)), since for each mode the model should be called once. Therefore the reduction step must reduce the number of modes as much as possible, but at the same time it must keep the important information in the covariance matrix. The advantage of this algorithm is that the original model operator is

used for the calculations and that the characteristic of a positive definite is preserved. This

algorithm is more robust but needs a large number of modes to be efficient.

(54)

44 4.3 Eulerian chemistry transport model EUROS

4.3 Eulerian chemistry transport model EUROS

The Eulerian atmospheric chemistry transport model EUROS [EURopean Operational Smog, Van Rheineck Leyssius et al. (1990); Van Pul et al. (1996); Matthijsen et al. (2001)], was devel-oped at RIVM (National Institute of Public Health and the Environment, The Netherlands). The model contains parameterizations of the various chemical, dynamical and radiative processes in the atmosphere, as well as information about the emissions of ozone precursors, nitrogen oxides, VOCs, carbon monoxide and methane. The model can be used to examine the time and

spatial behavior of , , , and in the lower troposphere over Europe.

4.3.1 Domain and processes

The model area extends over a large part of Europe and uses a shifted pole coordinate system,

i.e., the real North Pole is at northern latitude and the equator has been shifted to 600

northern latitude in the new coordinate system. In these shifted pole coordinates the model

domain is . The grid consists of grid cells with a

longitude-latitude resolution of (about 60 x 60 km). The vertical stratification of

Error subspace filtering for atmospheric chemistry data assimilation modeling

DATA ASSIMILATION MODELING

Proefschrift

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op maandag 6 november 2006 om 12.30 uur

door

Remus Gabriel HANEA

Contents

Chapter 1

Introduction

1.1

Air pollution

1.2

Ozone: "

"

1.3

Data assimilation

1.4

Overview

Chapter 2

Overview of the low-rank filter

algorithms

2.1

Introduction

2.2

State space model

2.3

Extended Kalman filter for nonlinear systems

2.3.1

Low-rank extended Kalman filters

2.3.2

Ensemble Kalman filters

Chapter 3

Convergence and divergence in the

Kalman filter approach

1

3.1

Introduction

3.2

The theory behind

3.3

Illustrative examples

3.4

Impact of suboptimal filter algorithms

3.4.1

Underestimation of the error covariance matrix

3.5

Illustrative examples (twin experiment)

3.5.1

The model

3.5.2

Uncertainties in emissions

3.6

Conclusions

Chapter 4

Kalman filter with a large scale

application

1

4.1

Introduction

4.2

Kalman filtering

4.2.1

State space model

4.2.2

Extended Kalman filter for nonlinear systems

4.2.3

Ensemble Kalman filter

4.2.4

Reduced Rank Square Root Kalman filter

4.3

Eulerian chemistry transport model EUROS

4.3.1

Domain and processes