Abstract

Temporally distributed deterministic and stochastic excitation of the tangent linear forecast system governing forecast error growth and the tangent linear observer system governing assimilation error growth is examined. The method used is to determine the optimal set of distributed deterministic and stochastic forcings of the forecast and observer systems over a chosen time interval. Distributed forcing of an unstable system addresses the effect of model error on forecast error in the presumably unstable forecast error system. Distributed forcing of a stable system addresses the effect on the assimilation of model error in the presumably stable data assimilation system viewed as a stable observer. In this study, model error refers both to extrinsic physical error forcing, such as that which arises from unresolved cumulus activity, and to intrinsic error sources arising from imperfections in the numerical model and in the physical parameterizations.

1. Introduction

A comprehensive understanding of the dynamics of linear systems is fundamental to physical science and, while linear systems arise in a wide variety of contexts, in this work we have in mind the particular examples of the dynamics of perturbations to the tangent linear forecast and data assimilation error systems. Among important problems in dynamical system theory is understanding the response of a linear system to forcing distributed in time. In the context of linear forecast error systems, temporally distributed forcing is used to account for the effect of model error on forecast predictability; while in the case of a data assimilation system it addresses the effect of model error on state identification for model initialization.

Because the forecast system is nonnormal as well as time dependent, its study requires the methods of generalized stability theory (GST; Farrell and Ioannou 1996a,b). These methods have heretofore been applied inter alia to the problem of cyclogenesis (Farrell 1982, 1989, the problem of transition to turbulence of boundary layer and laboratory channel flows (Farrell 1988a; Butler and Farrell 1992; Reddy and Henningson 1993), the maintenance of the turbulent state in laboratory shear flows (Farrell and Ioannou 1993b; Bamieh and Dahleh 2001; Jovanović and Bamieh 2001), the control of boundary layer and channel shear flow turbulence (Bewley and Liu 1998; Farrell and Ioannou 1998; Kim 2003; Högberg et al. 2003a,b), and in numerical weather prediction to determine the impact of uncertainties in the initial state on the forecast (Farrell 1988b; Lacarra and Talagrand 1988; Farrell 1990; Molteni and Palmer 1993; Buizza and Palmer 1995). However, as the initialization of forecast models is improved with the advent of new data sources and the introduction of variational assimilation methods, the medium range forecast will become increasingly affected by uncertainties resulting from incomplete physical parameterizations and numerical approximations in the forecast model, and by the necessarily misrepresented subgrid-scale chaotic processes such as cumulus convection that act as temporally stochastic distributed forcing of the forecast error system. These influences, referred to collectively as model error, conventionally appear as an external forcing in the forecast error system. Improving the understanding of model error and specifically identifying forcings that lead to the greatest forecast errors is centrally important in predictability studies. In analogy with the optimal perturbations that lead to the greatest forecast error in the case of initial condition error, these continuous error sources will be called optimal distributed forcings.

Two distinct dynamical representations for model error will be considered. The first is deterministic distributed forcing and the second is stochastic distributed forcing. The underlying theory for the deterministic case is based on analysis of the dynamical system as a mapping from the space of input forcings to the space of output states (Dullerud and Paganini 2000). This theory enabled Glover (1984) to obtain the optimal truncation of a dynamical system by using a singular value decomposition of the map between past forcings and future responses of the system (cf. Zhou and Doyle 1998). The singular values of this map, called Hankel singular values, have recently been used to truncate the linearized fluid equations and construct a reduced-order Kalman filter for use in meteorological assimilation (Farrell and Ioannou 2001a,b).

The link between distributed forcing and the growth of forecast errors is direct, but of equal interest is the role of distributed forcing in limiting the accuracy with which a system state can be identified. In order to address this problem we consider the data assimilation system viewed as a Luenberger observer (Luenberger 1971). Such an observer is necessarily asymptotically stable in contrast to the forward forecast problem, which is generally asymptotically unstable.

Recently D’Andrea and Vautard (2000) obtained approximate optimal temporally distributed deterministic forcings of the forecast error system (which they refer to as forcing singular vectors) and Barkmeijer et al. (2003) obtained the optimal temporally distributed deterministic forcing of the forecast error system over fixed spatial structures. In this paper, we extend the study of distributed forcings of the forecast error system by finding efficient methods for obtaining the exact solution of the optimal temporally distributed deterministic and stochastic forcing of the forecast and assimilation systems.

First we introduce the norms used to measure the forcing and the error state, and review the method used to obtain the optimal deterministic forcing that produces the greatest error state norm at a given time. The optimal problem for stochastic forcing and for impulsive forcing is then solved. Finally, examples are shown for optimal deterministic, stochastic, and impulsive forcing in forecast systems and observer systems.

2. Optimal deterministic forcing of linear systems

The problem of obtaining the forcing of a time-varying system over the time interval [0, T] that leads to the greatest state norm at time T has been solved in the context of control theory (see, i.e., Dullerud and Paganini 2000). Because of the novelty and power of the methods used to solve this problem, and their relation to methods used for obtaining optimal truncation of meteorological systems (cf. Farrell and Ioannou 2001a), it is useful to review this solution.

We seek the deterministic forcing f(t) of unit norm on t ∈ [0, T] producing the greatest state norm at time T that is, it maximizes the square norm of the state ||ψ(T)||2, assuming the state is initially zero, ψ(0) = 0, and that ψ obeys the linear equation:

 
formula

The forcing f(t) is taken to be in ℒ2([0, T]), that is, a square integrable vector function on the interval [0, T] with an inner product:

 
formula

The square norm in this space is

 
formula

The state ψ is in CN, that is, an N dimensional complex vector. We choose the Euclidean inner product for the state vectors:

 
formula

where H denotes the hermitian transpose. The square of the norm in this space is

 
formula

We have employed the subscript ℒ2 to denote the inner product of square integrable functions in order to distinguish it from the inner product of vectors.

With these norms we seek the forcing, f, maximizing the forcing normalized final state at time T:

 
formula

The state ψ at time T is obtained for any forcing, f, by the linear map ϕ:

 
formula

In (7) Φ(T, t) is the propagator associated with the linear dynamical system (1).

In order to obtain the optimal distributed forcing, some auxiliary results are needed. First the adjoint of the mapping ϕ is required. The adjoint of ϕ is denoted ϕ* and it maps a state ψ1(T) to a forcing g(t) defined over the interval [0, T], with the property that

 
formula

This adjoint map ϕ* is (cf. Farrell and Ioannou 2001a):

 
formula

Composition of the mapping ϕ and its adjoint ϕ* maps system states to system states. This is a rank N mapping that can be seen using (7) and (9) to have matrix representation:

 
formula

which is also the finite-time state covariance matrix under the assumption of temporally white noise forcing with unit covariance 𝗜. The covariance matrix 𝗖 in (10) is hermitian, and is assumed to be full rank. The covariance can also be obtained at time T by integrating from t = 0 to t = T the Lyapunov equation:

 
formula

with the initial condition 𝗖(0) = 0 (cf. Farrell and Ioannou 1996a).

Now we define the projection operator:

 
formula

Operator P maps forcings to forcings. This linear map is an orthogonal projector on the space of functions because it is hermitian and satisfies the property:

 
formula

Indeed,

 
formula

Also note the property:

 
formula

Indeed,

 
formula

Consider a forcing f, its projection Pf, and its orthogonal complement (IP)f. Then (15) establishes that forcing (IP)f makes no contribution to the state at ψ(T).

Assume that at time T, the system is in state ψ, and consider all the forcings f that could have produced this state, that is, all f with the property ϕf = ψ. That there exist forcings for which any chosen state ψ can be reached at time T requires that the system be controllable, but this is ensured by the full rank of 𝗖, which in the control literature is referred to as the controllability Grammian (Brockett 1970).

All the forcings, f, that produce state ψ have the same projection Pf. To show this, note that by assumption:

 
formula

and by (15) the projected forcing Pf also produces the state ψ so

 
formula

Multiplying (18) by ϕ*𝗖−1 gives

 
formula

But the left side of (19) is P2f = Pf. Hence all f that produce state ψ have the same projection:

 
formula

It is now easy to see that the forcing Pf is the minimal forcing producing the state ψ. To show this, we decompose the forcing into its orthogonal complements:

 
formula

Because

 
formula

we have

 
formula

and it follows that

 
formula

which proves that the forcing

 
formula

is the minimal forcing to produce state ψ. It is easy to verify that this forcing produces ψ:

 
formula

Using (24) the maximization (6) over functions, f, can be transformed to a maximization over states, ψ(T):

 
formula

This equation shows that the optimal over forcings (lhs) is equal to an optimal over states (rhs). Note that the form of the second optimization is reminiscent of the covariance or Mahalanobis metric often used in predictability analysis (Palmer et al. 1998) suggesting the interpretation of optimals weighted by the Mahalanobis metric as structures that are most easily forced.

Quotient (27) is maximized for unit forcing by the state:

 
formula

where λ1 is the maximum singular value of 𝗖 and υ1 is the corresponding singular vector of 𝗖 [υ1 is conventionally called the top empirical orthogonal function (EOF) of 𝗖]. The optimal unit forcing is obtained from (9) and (25):

 
formula

and the optimal value of the quotient (27) is clearly λ1. It is also clear from (29) that the optimal forcing structure at optimizing time, fopt(T), is identical to the optimal state υ1. If the system is asymptotically stable, λ1 remains finite as T → ∞, and if the system is autonomous as well as stable, λ1 asymptotes to a constant; while if the system is asymptotically unstable, λ1 grows exponentially as T → ∞. In the limit of small optimizing time 𝗖 = O(T) and therefore (27) approaches zero linearly in time, a result that can be traced to the quadratic weighting of the forcing by ||·||22, which has the effect of penalizing concentration of forcing in time.

If υ1 and υ2 are two distinct singular vectors of 𝗖, the corresponding forcings f1,2(t) = ϕ*𝗖−1υ1,2 are orthogonal because

 
formula

The optimal forcing and the associated state of the system can be obtained simultaneously by integrating the following coupled forward and adjoint systems backward over the finite interval from time t = T to the initial time t = 0:

 
formula

with ψ(T) = λ1υ1 and f(T) = υ1/λ1. The initial state ψopt(0) = 0 is recovered as a consistency check.

Consider forcings, f, with ||f||22 = 1. The states these produce will lie on an ellipsoid with axes of length λi in the directions of the orthonormal υi, where υi are the singular vectors of 𝗖 with corresponding singular values λi. The corresponding forcings

 
formula

which produce the extremal states ψi(T) = λiυi are an orthonormal finite dimensional set of functions that span the space of minimal forcings. These forcings induce a characterization of ψ(T) and its excitation by distributed forcing analogous to the characterization induced by the evolved optimals and the optimal initial conditions in the initial condition excitation problem. For this reason we extend the characterization of optimality to all the forcings that render the Rayleigh quotient (6) extremal. A distinction between optimal distributed forcings and optimal initial conditions is that the set of optimal distributed forcings, fi, is the minimal spanning set of forcings that produce nonzero states at t = T, but the set of optimal forcings, fi, do not span the whole space ℒ2([0, T]). However, the optimal forcings do split the whole of ℒ2([0, T]) into two orthogonal components: the range of P, P2([0, T]), comprises linear combinations of the orthogonal set of fi which span the forcings that contribute to the state at T, and the orthogonal space of functions in the kernel of P, (IP)ℒ2([0, T]), which do not contribute to the state at T.

An arbitrary forcing f may be projected on its minimal component Pf with the property that only this component of the forcing produces a nonzero state at the optimizing time T.

Specifically, an arbitrary state ψ(T) can be expanded in the basis of the orthonormal eigenvectors, υi, of the finite-time covariance matrix 𝗖(T):

 
formula

where

 
formula

Then the minimal forcing,

 
formula

is the only relevant forcing for producing this state. Further, because λi generally fall rapidly with i, most states are difficult to reach; for instance by (35), states ψ(T) that project on υi with small values of λi require large-amplitude forcing. It follows that calculation of the first few eigenvectors, υi, of the finite-state covariance 𝗖(T) and the optimal forcings, fi, that produce them, suffices to characterize the potential influence of model error on the forecast and assimilation systems. In this way the decomposition of forcings into contributing and noncontributing spaces can be made independent of the state dimension N, when N is large, by taking the spanning set for the contributing space to be restricted to the forcings with substantial λi.

3. Optimal constant spatial structure for stochastic forcing of a linear system

In the previous section we obtained the optimal deterministic forcing for producing a specified state and the maximum state norm at a given time. However, situations arise in which the forcing is not known precisely but can be characterized statistically. A common characteristic of statistically specified, physically derived forcing is that the forcing is stochastic with a short correlation time when compared to the characteristic time scales of the linear operator. Given such a stochastic time dependence we seek the constant spatial forcing structure producing the greatest expected state square norm at a chosen time. This situation is formulated using the following stochastic equation:

 
formula

where f is a column vector specifying the spatial structure and η(t) a temporally white noise scalar function of unit variance satisfying 〈η(t1)η(t2)〉 = δ(t2t1).

The optimal spatial structure, fopt, of the stochastic forcing is to be found under the assumption that its temporal behavior is a white noise process, η(t). This contrasts with the deterministic optimal forcing problem studied above, in which the optimal response was found over both the spatial and temporal structure of the forcing.

The problem of obtaining the fixed spatial structure of temporally white forcing that excites the greatest expected square state norm in the asymptotic limit t → ∞ has been addressed previously (Farrell and Ioannou 1993a,b, 1996a; Kleeman and Moore 1997). This structure is the stochastic optimal. The problem of obtaining the analogous finite-time stochastic optimal under both certain and uncertain dynamics was studied by Farrell and Ioannou (2002a,b). We will now briefly review how the stochastic optimal is determined.

The ensemble expected square norm of the state at time T is

 
formula
 
formula

where the matrix 𝗕(T) is

 
formula

The matrix 𝗕(T) can be obtained by integrating from t = 0 to t = T the equation,

 
formula

with the initial condition 𝗕(0) = 0.

It is immediate from (38) that the singular vectors of the hermitian matrix 𝗕 are the finite-time stochastic optimals. The finite-time stochastic optimal associated with the largest singular value of 𝗕 is the unit norm spatial forcing structure that, when excited white in time over the interval [0, T], produces the greatest expected square norm of the state at time T. This unit norm spatial structure results in unit variance forcing in (36). This largest expected square state response, when normalized by the forcing covariance, is

 
formula

where μ1 is the largest singular value of 𝗕(T) and fopt the corresponding singular vector. The forcing normalizations in (41) are not the same as in the deterministic case. Deterministic forcings are normalized by the forcing’s integrated square amplitude as in (3). Such a normalization is not appropriate for white noise forcing because white noise realizations are unbounded in square amplitude. The normalization in (41) assures a unit variance forcing in (36), which is appropriate for white noise forcing.

4. Optimal impulsive forcing on an interval

Consider an impulsive forcing applied at a point in time, t0:

 
formula

where f is the spatial structure of the forcing. It is assumed, as previously, that the state is initially zero, ψ(0) = 0, and that the forcing is applied at a point in the interval so that, 0 ≤ t0T, where T is a final time.

We seek the time t0 and the forcing spatial structure f of unit norm, corresponding to unit impulse, resulting in the greatest state norm at T, that is, we seek to maximize the quotient:

 
formula

The state at time T is

 
formula

Consequently, we seek the spatial structure that maximizes

 
formula

The maximum value of Rp is the square of the largest singular value σ21, of the propagator Φ(T, t0), which is also the square of the norm of Φ(T, t0), denoted ||Φ(T, t0)||22. The forcing applied at t0 producing the greatest state norm at T is the top singular vector of the propagator Φ(T, t0), which is the deterministic global optimal initial condition corresponding to optimizing time T − t0. Therefore, the optimal point forcing has the structure of the global optimal initial condition over t ∈ [0, T] and is applied at the global optimal time in the past. If we denote by tgo the time t for which ||Φ(T, t)||2 is maximum, that is, the time at which the global optimal growth is achieved, then the forcing is applied at t0 = tgo. It is clear that the optimal point forcing is identical to the optimal initial condition at t0.

Given an arbitrary final state ψ(T) the optimal impulsive forcing to produce ψ(T) is obtained by finding

 
formula

in which Φ(t, T) = Φ−1(T, t), and taking t0 as the time of this minimum and f = Φ(t, T)ψ(T)δ(tt0) as the impulsive forcing that optimally excites state ψ(T).

A related problem is finding the minimum constant amplitude time-varying structure over an interval that produces a given state. The solution is to force the normalized function obtained from the backward propagator:

 
formula

5. Examples

a. A scalar system example

Consider the one-dimensional and necessarily normal system:

 
formula

From (28) and (29) and the fact that a one-dimensional system is self-adjoint, it follows that the deterministic unit norm forcing maximizing the state at optimizing time T is

 
formula

and the corresponding maximum state at optimizing time T is

 
formula

The optimal forcing for this example is shown in Fig. 1a for a representative stable case (α = −1) and in Fig. 1b for a similar unstable case (α = 1); both for T = 5. Note that the forcing in the stable case is concentrated near the optimizing time at the end of the interval, because forcing at earlier times would excite the decaying state and, as a result, be suboptimal. In the unstable case the forcing is concentrated at the start of the interval to take advantage of the growth. This is a general property of optimal forcings. The square norm of the state at the optimizing time T as a function of T is shown for the stable case in Fig. 1c (case labeled deterministic). For small T the square norm of the state is |ψopt(T)|2T but for large T it asymptotes to 1/2α. The optimal forcing at optimizing time T has norm fopt(T) = 1/λ while the optimal state norm is ψopt = λ with an optimal finite-time variance:

 
formula
Fig. 1.

For the one-dimensional system (48): (a) unit norm optimal forcing as a function of time for the stable case with α = −1 and T = 5, (b) unit norm optimal forcing for the unstable case with α = 1 and T = 5, (c) square of the final state for the optimal deterministic, stochastic, and impulsive forcings at optimizing time T as a function of the optimizing time for the stable case with α = −1. The optimal deterministic and the optimal stochastic responses coincide. (d) A decomposition of f(t) = 1 into the component contributing, Pf, and the component not contributing, fPf, to the state at t = 5; also for α = −1

Fig. 1.

For the one-dimensional system (48): (a) unit norm optimal forcing as a function of time for the stable case with α = −1 and T = 5, (b) unit norm optimal forcing for the unstable case with α = 1 and T = 5, (c) square of the final state for the optimal deterministic, stochastic, and impulsive forcings at optimizing time T as a function of the optimizing time for the stable case with α = −1. The optimal deterministic and the optimal stochastic responses coincide. (d) A decomposition of f(t) = 1 into the component contributing, Pf, and the component not contributing, fPf, to the state at t = 5; also for α = −1

As in this example, it is generally the case that in stable normal systems the optimal forcing extends over a time interval of the order of the decay time of the system while in nonnormal systems the optimal forcing assumes the temporal extent of the system’s global optimal as will be shown. The expected square norm of the state at T for unit covariance stochastic forcing is easily seen to be

 
formula

which is identical to the optimal deterministic state square norm. This is shown in Fig. 1 for a = −1. Also shown in Fig. 1 is the norm of the state for optimal impulsive forcing, which is identical to the optimal obtained from any unit initial condition on the interval. Because in the case of a stable normal system the amplitude of the states decays monotonically after they are excited, the optimal point for impulsive forcing is at the optimizing time T itself and consequently the optimal response is equal to 1 for all choices of T.

Assume that in this simple scalar system we are given a forcing f(t). To find the part of f(t) that influences the state at T, we need only to calculate the projection of f(t) on the space of optimal forcing that according to (12) is

 
formula

The part of f(t) that does not influence the final state is fPf. For example if the forcing were constant, f(t) = 1, only the projection,

 
formula

produces a contribution to the state at T, while fPf produces no contribution. This decomposition of f(t) = 1 into contributing and noncontributing components is shown in Fig. 1.

b. The Reynolds system example

Consider the two-dimensional highly nonnormal Reynolds system that models streak formation in boundary layers:

 
formula

with stable matrix 𝗔:

 
formula

The square norm of the state at time T produced by optimal deterministic forcing over the interval t ∈ [0, T] is shown as a function of T in Fig. 2. The response is much higher than that found for a normal system with the same eigenvalues because of the nonnormality of the system [e.g., (56)]. The variance equation for the above system is

 
formula

The first term ψH(𝗔H + 𝗔)ψ represents the variance tendency from the background mediated by the nonnormality of 𝗔 while (fψH + ψfH) is the variance tendency directly accumulated from the forcing. These two terms are plotted in Fig. 3c as a function of time for the forcing that leads to maximum state norm at time T = 5 in system (55). Note that the variance obtained directly from forcing is far less than the variance obtained from nonnormal processes.

Fig. 2.

For the Reynolds system example: the optimal square norm of the state ||ψopt(T)||2 at optimizing time for deterministic, stochastic, and impulsive forcing as a function of optimizing time T. For comparison, the optimal growth for unit norm initial conditions, ||eAT||2, is also shown. Note that the optimal stochastic response coincides with the optimal deterministic response in this two-dimensional example

Fig. 2.

For the Reynolds system example: the optimal square norm of the state ||ψopt(T)||2 at optimizing time for deterministic, stochastic, and impulsive forcing as a function of optimizing time T. For comparison, the optimal growth for unit norm initial conditions, ||eAT||2, is also shown. Note that the optimal stochastic response coincides with the optimal deterministic response in this two-dimensional example

Fig. 3.

For the Reynolds system example: (a) time development of the components of the optimal forcing f = (f1, f2) that optimize the state energy at T = 5; (b) corresponding evolution of the square norm of the state; (c) variance tendency due to nonnormal processes ψHopt(t) (𝗔 + 𝗔H)ψopt(t) and to forcing fHopt(t)ψopt(t) + ψHopt(t)fopt(t) as a function of time. The variance tendency due to nonnormal processes exceeds variance tendency due to forcing

Fig. 3.

For the Reynolds system example: (a) time development of the components of the optimal forcing f = (f1, f2) that optimize the state energy at T = 5; (b) corresponding evolution of the square norm of the state; (c) variance tendency due to nonnormal processes ψHopt(t) (𝗔 + 𝗔H)ψopt(t) and to forcing fHopt(t)ψopt(t) + ψHopt(t)fopt(t) as a function of time. The variance tendency due to nonnormal processes exceeds variance tendency due to forcing

Because this is a two-dimensional example, the state at t = T produced by the optimal stochastic and deterministic forcing are identical.1 The optimal norm of the state produced by impulsive forcing and, for comparison, the optimal initial value growth ||e𝗔T||2 are shown as a function of T in Fig. 2. Note that the norm of the state produced by optimal impulsive forcing follows the optimal growth curve until T reaches the global optimal time tgo (which occurs for this case at tgo = 0.67), then for T > tgo it maintains this constant value. The optimal impulsive forcing is applied at time t = Ttgo and has the structure of the global optimal initial condition.

The structure of the forcing as a function of time is shown in Fig. 3a and the corresponding evolution of the state norm is shown in Fig. 3b. The forcing is seen to peak tgo units of times prior to the optimizing time T.

6. Modeling data assimilation as a stable observer system

Consider assimilating data taken from the exact solution of a nonlinear system, xt. The forecast error ef = xfxt obeys the following equation:

 
formula

in which model error is represented by f(t) and 𝗔 is the tangent linear system about xt. Forecast error typically grows exponentially when averaged over a sufficiently long time due to the Lyapunov instability that inherently results from the time dependence of 𝗔 (Arnold 1998; Farrell and Ioannou 1999). However, a Lyapunov stable system related to (58) can be constructed to model the process of state estimation for forecast initialization, a process referred to as data assimilation. This can be done because a data assimilation system is equivalent to a system in which the matrix 𝗔 in (58) is modified by the introduction of observations so that the solution to the modified system, xa, converges to the system state, xt, in a statistical sense (Luenberger 1971). The effect of model and observational error on data assimilation can be analyzed quite generally using this representation of a data assimilation system as a Luenberger observer.

Assume that observations, yob of the state are taken at various locations and introduced into the model. Let the observations be related to the true state of the atmosphere, xt, as

 
formula

where 𝗛 is an observation matrix, 𝗥 is the observational error covariance, and w is a vector of white noise processes. Given an observational platform, yob is available while the truth xt is not available.

Assimilate these observations to obtain an analysis, xa, with analysis error ea = xaxt satisfying the Luenberger observer system:

 
formula

The above equation defines the assimilation system model. Ideally, the gain, 𝗞, is chosen to minimize the analysis error variance trace (〈eaeHa〉).

Unlike the forecast error system, a Luenberger observer system is asymptotically stable. Any gain, 𝗞, that stablizes 𝗔 − 𝗞𝗛 results in an observer with bounded error, this error being forced by a combination of model error f(t) and observational error 𝗥 (cf. 60). It is immediately evident that in the absence of model and observation error the analysis will converge to truth. Good gains, 𝗞, do not just stabilize the operator but simultaneously reduce the nonnormality of the tangent linear operator so that the minimum of trace (〈eaeHa〉) is maintained by the combination of observational and model error. The optimal gain under the assumption of Gaussian error statistics is the Kalman gain (Kalman 1964; Zhou and Doyle 1998). A suboptimal observer, such as might be associated with an operational forecast system, can be modeled as a stabilizing 𝗞 not equal to the Kalman gain. Just as analysis of the tangent linear forecast system reveals how initialization error and distributed forcing contribute to forecast error, so does analysis of the observer system reveal how model and observation error contribute to analysis error.

We are now equipped to study the effect of distributed error on data assimilation and will begin by seeking the deterministic and stochastic forcings leading to the greatest deviations of the estimate xa from the true state xt as measured by the norm of e.

7. The Eady assimilation and forecast system model example

Consider as an example error system the Eady model (Eady 1949) with a zonal wavenumber chosen so that the model is unstable. Introduce continuous observation of the error streamfunction at the channel center and construct a stable Kalman observer under the assumption of zero model and unit observation error. Under these assumptions, a single continuous observation at the channel center suffices to make this optimal observer asymptotically stable.

The nondimensional linearized equation for streamfunction perturbations ψ in the Eady model is

 
formula

in which the perturbation is assumed to be of form ψ(z, t)eikx, where k is the zonal wavenumber. The perturbation potential vorticity is D2ψ, with D2 ≡ ∂2/∂z2k2 and the perturbation potential vorticity damping rate is r. This model is Boussinesq with constant stratification on an f plane and has periodic boundary conditions in the zonal, x, direction and a solid lid at height z = H. The zonal flow is U(z) = z. Horizontal scales are nondimensionalized by L = 1000 km; vertical scales by H = fL/N = 10 km; velocity by U0 = 30 m s; and time by T = L/U0, so that a time unit is approximately 9.25 h. Midlatitude values are chosen for the Brunt–Väisälä frequency, N = 10−2 s−1, and Coriolis parameter, f = 10−4 s−1. The coefficient of Rayleigh damping has value r = 0.15 corresponding to a decay rate of 1/2.5 day−1.

Conservation of potential temperature at the ground, z = 0, and tropopause, z = 1, provides the following boundary conditions:

 
formula
 
formula

The model is discretized using N = 21 points in the vertical and expressed in the matrix form:

 
formula

In (64) generalized velocity coordinates ϕ = 𝗠1/2ψ were used so that the Euclidian square norm of the system state corresponds to error energy. The error energy is ϕHϕ = ψH𝗠ψ, and for a vertical discretization increment δz the metric 𝗠 is

 
formula

The assimilation system in generalized velocity variables is governed by the operator 𝗔a = 𝗔 − 𝗞𝗛, where 𝗞 is the Kalman gain for observation at the channel center given by 𝗞 = 𝗫𝗛H𝗥−1, where 𝗫 satisfies the continuous-time algebraic Ricatti equation under the perfect model assumption (cf. Zhou and Doyle 1998):

 
formula

and 𝗛 = 𝗣𝗠−1/2 is the observation matrix, in which, for our example, 𝗣 is a row vector with one nonzero element corresponding to the grid point at the channel center, and 𝗥 is the observation error covariance matrix (which is a scalar in our example). In contrast, the unstable forecast system is governed by the operator 𝗔. We seek in both the stable assimilation and the unstable forecast system the optimal deterministic and stochastic forcing leading to the greatest error at time T.

a. The observer-stabilized Eady model

Let 𝗔a be the matrix of the observer-stabilized Eady model for zonal wavenumber k = 1.5. While the unobserved Eady model at this wavenumber is unstable with dimensional growth rate 1/2.5 day−1, the observer-stabilized model has a minimum dimensional decay rate 1/2.6 day−1. The square norm of the state at time T produced by optimal distributed forcing over the interval t ∈ [0, T] is shown as a function of T in Fig. 4. The square norm of the state at T produced by each of the 21 orthonormal optimal forcings is shown in Fig. 5. Two EOFs and their associated forcing structures account for most of the perturbation dynamics suggesting the generally valid conclusion that a few forcings account for most of the error and the majority of the forcings make relatively little contribution.

Fig. 4.

For the observer-stabilized Eady model: the square norm of the optimal state at optimizing time T for unit norm deterministic distributed forcing; the optimal expected square norm of the state for statistically steady stochastic forcing; and the optimal square norm of the state for impulsive forcing. The zonal wavenumber is k = 1.5 and the Rayleigh damping coefficient is r = 0.15. The Eady model is stabilized by the asymptotic Kalman gain, obtained under the perfect model assumption and with unit observation error variance. The unobserved system is unstable with dimensional growth rate 1/2.3 day−1. The decay rate of perturbations in the observer system is 1/2.6 day−1

Fig. 4.

For the observer-stabilized Eady model: the square norm of the optimal state at optimizing time T for unit norm deterministic distributed forcing; the optimal expected square norm of the state for statistically steady stochastic forcing; and the optimal square norm of the state for impulsive forcing. The zonal wavenumber is k = 1.5 and the Rayleigh damping coefficient is r = 0.15. The Eady model is stabilized by the asymptotic Kalman gain, obtained under the perfect model assumption and with unit observation error variance. The unobserved system is unstable with dimensional growth rate 1/2.3 day−1. The decay rate of perturbations in the observer system is 1/2.6 day−1

Fig. 5.

Optimal deterministic square state norm at T = 40 in the observer-stabilized Eady model accounted for by each of the 21 orthogonal forcing structures. There are two dominant responses to deterministic forcing, suggesting the generally valid result that a few forcings account for most assimilation error. The parameters are as described in Fig. 4 

Fig. 5.

Optimal deterministic square state norm at T = 40 in the observer-stabilized Eady model accounted for by each of the 21 orthogonal forcing structures. There are two dominant responses to deterministic forcing, suggesting the generally valid result that a few forcings account for most assimilation error. The parameters are as described in Fig. 4 

The norm of the distributed forcing that maximizes the state norm at T = 40 is shown as a function of time in Fig. 6a along with the optimal growth ||e𝗔a(Tt)||. The norm of the forcing follows the optimal growth curve and is maximized at T − tgo, where tgo is the time for which the maximum of the norm of the propagator e𝗔at occurs (in this case the maximum occurs at tgo = 7). The corresponding growth of the error state energy, the energy tendency by nonnormal processes, and the energy tendency produced by forcing are shown in Fig. 6b.

Fig. 6.

For the observer-stabilized Eady model: (a) time development of the norm of the deterministic forcing that optimizes the energy of the state at T = 40. The norm of the retarded propagator ||e𝗔a(Tt)|| giving the optimal growth attainable from an initial state is also shown. Note that the optimal deterministic forcing is concentrated in the interval of the supreme of the maximum norm of the retarded propagator. (b) The corresponding time development of the norm of the state ||ψopt(t)||, the variance tendency due to nonnormal processes ψHopt(t)(𝗔 + 𝗔H)ψopt(t), and the variance tendency due to forcing fHopt(t)ψopt(t) + ψHopt(t)fopt(t). The forcing is normalized so that ||fopt(T)||2 = 1. The parameters are the same as in Fig. 4 

Fig. 6.

For the observer-stabilized Eady model: (a) time development of the norm of the deterministic forcing that optimizes the energy of the state at T = 40. The norm of the retarded propagator ||e𝗔a(Tt)|| giving the optimal growth attainable from an initial state is also shown. Note that the optimal deterministic forcing is concentrated in the interval of the supreme of the maximum norm of the retarded propagator. (b) The corresponding time development of the norm of the state ||ψopt(t)||, the variance tendency due to nonnormal processes ψHopt(t)(𝗔 + 𝗔H)ψopt(t), and the variance tendency due to forcing fHopt(t)ψopt(t) + ψHopt(t)fopt(t). The forcing is normalized so that ||fopt(T)||2 = 1. The parameters are the same as in Fig. 4 

The evolution in structure of the first deterministic optimal is shown in the left panels of Fig. 7. The corresponding evolution in structure of the state is shown in the right panels. At optimizing time both the state and the optimal forcing assume the structure of the first eigenvector of the finite-time covariance matrix at time T, 𝗖(T).

Fig. 7.

For the observer-stabilized Eady model: (left) temporal evolution of the structure of the deterministic optimal forcing Re [f(z, t)eikx] for T = 40. The optimal forcing is shown from the bottom to the top at times t = 20, 30, 35, 40. (right) The corresponding structure of the optimal state Re [ψ(z, t)eikx]. The parameters are as described in Fig. 4 

Fig. 7.

For the observer-stabilized Eady model: (left) temporal evolution of the structure of the deterministic optimal forcing Re [f(z, t)eikx] for T = 40. The optimal forcing is shown from the bottom to the top at times t = 20, 30, 35, 40. (right) The corresponding structure of the optimal state Re [ψ(z, t)eikx]. The parameters are as described in Fig. 4 

The error norm produced by optimal stochastic forcing, and the response to optimal impulsive forcing are shown as a function of T in Fig. 4. The norm of the error due to impulsive forcing follows the deterministic initial state optimal growth curve reaching at time tgo = 7 the global optimal norm after which it assumes the constant value of the global optimal growth. For optimizing times T > tgo the optimal impulsive forcing occurs at time t = Ttgo and has the structure of the global optimal initial condition. This structure and that of the stochastic optimal are shown in Fig. 8. If the spatial structure of the white noise forcing is taken to be that of the stochastic optimal shown in Fig. 8a then the resulting state covariance at optimizing time is primarily composed of the first three eigenvectors of 𝗖(T) shown in Fig. 9.

Fig. 8.

For the observer-stabilized Eady model: (a) structure of the stochastic optimal that produces the greatest expected square error at T = 40, (b) structure of the optimal impulsive forcing that produces the greatest expected square error at T = 40. The optimal impulsive forcing is concentrated 7 units of time before the optimization time and has the structure of the optimal initial condition associated with the global optimal growth. The parameters are as described in Fig. 4 

Fig. 8.

For the observer-stabilized Eady model: (a) structure of the stochastic optimal that produces the greatest expected square error at T = 40, (b) structure of the optimal impulsive forcing that produces the greatest expected square error at T = 40. The optimal impulsive forcing is concentrated 7 units of time before the optimization time and has the structure of the optimal initial condition associated with the global optimal growth. The parameters are as described in Fig. 4 

Fig. 9.

Structure of the (left to right) first three eigenvectors of 𝗖(T) (EOFs) at T = 40 maintained by the stochastic optimal in the observer-stabilized Eady model. The parameters are as described in Fig. 4 

Fig. 9.

Structure of the (left to right) first three eigenvectors of 𝗖(T) (EOFs) at T = 40 maintained by the stochastic optimal in the observer-stabilized Eady model. The parameters are as described in Fig. 4 

b. The unstable Eady forecast model

The error in the unstable Eady model at k = 1.5 grows exponentially. The error norm for optimal distributed, stochastic, and impulsive forcing are shown in Fig. 10. The optimal stochastic forcing and the optimal deterministic forcing produce almost identical error at time T because the system is unstable and excitation of a single structure, the adjoint of the mode, is responsible for nearly all the error at time T. The time development of the norm of the optimal deterministic forcing for optimizing time T = 20 is shown in Fig. 11. The forcing is concentrated at the beginning of the interval in order to take advantage of the instability. The structure of the forcing and the error at various times are shown in Fig. 12. At the optimizing time the structure both of the error and the forcing is that of the fastest growing mode, at earlier times the forcing has the structure of the adjoint of this mode. The stochastic optimal and the optimal impulsive forcing (which consists of forcing at the initial time of and with the structure of the optimal initial condition that maximizes perturbation growth at T = 20), although formally different, in this strongly unstable example have both assumed the structure of the adjoint of the most unstable mode of the model, which is shown in the second panel from the top on the left-hand side of Fig. 12.

Fig. 10.

For the unstable Eady forecast error model: the square norm of the state, ||ψopt(T)||2, at optimizing time, T, as a function of optimizing time under deterministic forcing of unit norm, which, for this strongly unstable error model, coincides with the optimal expected square error norm under statistically stationary stochastic forcing. For comparison, the optimal square error norm under optimal impulsive forcing is also shown. The zonal wavenumber is k = 1.5 and the Rayleigh damping coefficient is r = 0.15. The forecast system is unstable with a growth rate 1/2.5 day−1

Fig. 10.

For the unstable Eady forecast error model: the square norm of the state, ||ψopt(T)||2, at optimizing time, T, as a function of optimizing time under deterministic forcing of unit norm, which, for this strongly unstable error model, coincides with the optimal expected square error norm under statistically stationary stochastic forcing. For comparison, the optimal square error norm under optimal impulsive forcing is also shown. The zonal wavenumber is k = 1.5 and the Rayleigh damping coefficient is r = 0.15. The forecast system is unstable with a growth rate 1/2.5 day−1

Fig. 11.

For the unstable Eady forecast model example: (a) norm of the deterministic forcing that optimizes the energy of the state at T = 20. (b) Corresponding norm of the state |ψopt(t)|. The parameters are the same as described in Fig. 10 

Fig. 11.

For the unstable Eady forecast model example: (a) norm of the deterministic forcing that optimizes the energy of the state at T = 20. (b) Corresponding norm of the state |ψopt(t)|. The parameters are the same as described in Fig. 10 

Fig. 12.

For the unstable Eady forecast model: (left) structure of the deterministic optimal forcing Re [f(z, t)eikx] that optimizes the state norm at T = 20. The optimal forcing is shown at times t = 2, 10, 15, 20. (right) the corresponding structure of the state Re [ψ(z, t)eikx]. The parameters are as described in Fig. 10 

Fig. 12.

For the unstable Eady forecast model: (left) structure of the deterministic optimal forcing Re [f(z, t)eikx] that optimizes the state norm at T = 20. The optimal forcing is shown at times t = 2, 10, 15, 20. (right) the corresponding structure of the state Re [ψ(z, t)eikx]. The parameters are as described in Fig. 10 

8. Conclusions

In this work we have studied the response of a linear dynamical system to distributed deterministic and stochastic forcing. Specific atmospheric applications of this study are to the presumed unstable forward forecast error system and the necessarily stable data assimilation system that are each continually perturbed by disturbances that degrade the forecast and assimilation, respectively. These disturbances may be physical in origin and may also arise from numerical model deficiencies. The disturbances may be deterministic or stochastic. We analyzed the effect of temporally distributed forcing on stable and unstable error systems by obtaining a spanning set of optimal distributed forcings over a chosen interval making use of the adjoint system propagator applied to the finite-time covariance matrix eigenvectors. In the stable case modeling a data assimilation system, the optimal forcing is concentrated near the end of the interval of the optimization, while in the unstable forecast error system the optimal forcing is concentrated near the beginning of the interval. The optimals use both accumulation of forcing variance as well as variance transfer from the mean in their growth process, although transfer from the mean was dominant in the examples. The structure of the response to the optimal distributed forcing is the first eigenvector of the finite-time covariance 𝗖(T) = ∫T0Φ(T, t)ΦH(T, t)dt. In the case of stochastic forcing, the structure of the forcing producing the greatest response is the finite-time stochastic optimal, which is the first eigenvector of the stochastic optimal matrix 𝗕(T) = ∫T0ΦH(T, t)Φ(T, t)dt, and this forcing generally tends to have the structure of the optimal initial condition. Finally, in the case of point impulsive forcing, the greatest response is obtained by a forcing concentrated at the global optimizing time in the past and with the structure of the global optimal initial condition.

Most forcings make little or no contribution to forecast and assimilation error, and a small subspace of forcings produce most of the error that is realized primarily in a small subspace of error states. These forcings are the optimal forcings that project on the dominant eigenvectors of the finite-time covariance matrix.

Forecast error may arise primarily from error in the initial analysis, or it may be forced by error sources distributed over the forecast interval. Initial error influence on forecast is characterized by the set of optimal perturbations, while error distributed in time is characterized by optimal forcing functions. It remains to determine the relative importance of these error sources in practical forecast situations.

Consider macroscopic cyclone formation as distinct from perturbation error growth. The cyclone may arise from forcing distributed over a given interval or from an initial disturbance at the beginning of the interval. Improved understanding of the origin of cyclones depends, in part, on a better characterization of the relative role of initial conditions and distributed forcings in the cyclogenesis process.

Acknowledgments

This work was supported by NSF Grant ATM-0123389 and by ONR Grant N00014-99-1-0018.

REFERENCES

REFERENCES
Arnold
,
L.
,
1998
:
Random Dynamical Systems.
Springer-Verlag, 586 pp
.
Bamieh
,
B.
, and
M.
Dahleh
,
2001
:
Energy amplification in channel flows with stochastic excitation.
Phys. Fluids
,
13
,
3258
3269
.
Barkmeijer
,
J.
,
T.
Iversen
, and
T. N.
Palmer
,
2003
:
Forcing singular vectors and other sensitive model structures.
Quart. J. Roy. Meteor. Soc.
,
129
,
2401
2424
.
Bewley
,
T. R.
, and
S.
Liu
,
1998
:
Optimal and robust control and estimation of linear paths to transition.
J. Fluid Mech.
,
365
,
23
57
.
Brockett
,
R. W.
,
1970
:
Finite Dimensional Linear Systems.
John Wiley and Sons, 244 pp
.
Buizza
,
R.
, and
T. N.
Palmer
,
1995
:
The singular-vector structure of the atmospheric global circulation.
J. Atmos. Sci.
,
52
,
1434
1456
.
Butler
,
K. M.
, and
B. F.
Farrell
,
1992
:
Three dimensional optimal perturbations in viscous shear flow.
Phys. Fluids,
,
4A
,
1637
1650
.
D’Andrea
,
F.
, and
R.
Vautard
,
2000
:
Reducing systematic errors by empirically correcting model errors.
Tellus
,
52A
,
21
41
.
Dullerud
,
G. E.
, and
F.
Paganini
,
2000
:
A Course in Robust Control Theory: A Convex Approach.
Springer-Verlag, 440 pp
.
Eady
,
E. A.
,
1949
:
Long waves and cyclone waves.
Tellus
,
1
,
33
52
.
Farrell
,
B. F.
,
1982
:
The initial growth of disturbances in a baroclinic flow.
J. Atmos. Sci.
,
39
,
1663
1686
.
Farrell
,
B. F.
,
1988a
:
Optimal excitation of perturbations in viscous shear flow.
Phys. Fluids
,
31
,
2093
2102
.
Farrell
,
B. F.
,
1988b
:
Optimal excitation of neutral Rossby waves.
J. Atmos. Sci.
,
45
,
163
172
.
Farrell
,
B. F.
,
1989
:
Optimal excitation of baroclinic waves.
J. Atmos. Sci.
,
46
,
1193
1206
.
Farrell
,
B. F.
,
1990
:
Small error dynamics and the predictability of atmospheric flows.
J. Atmos. Sci.
,
47
,
2409
2416
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1993a
:
Stochastic dynamics of baroclinic waves.
J. Atmos. Sci.
,
50
,
4044
4057
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1993b
:
Stochastic forcing of the linearized Navier–Stokes equations.
Phys. Fluids,
,
5A
,
2600
2609
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1996a
:
Generalized stability. Part I: Autonomous operators.
J. Atmos. Sci.
,
53
,
2025
2041
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1996b
:
Generalized stability. Part II: Nonautonomous operators.
J. Atmos. Sci.
,
53
,
2041
2053
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1998
:
Turbulence suppression by active control.
Phys. Fluids
,
8
,
1257
1268
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
1999
:
Perturbation growth and structure in time-dependent flows.
J. Atmos. Sci.
,
56
,
3622
3639
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2001a
:
Accurate low-dimensional approximation of the linear dynamics of fluid flow.
J. Atmos. Sci.
,
58
,
2771
2789
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2001b
:
State estimation using a reduced-order Kalman filter.
J. Atmos. Sci.
,
58
,
3666
3680
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2002a
:
Perturbation growth and structure in uncertain flows. Part II.
J. Atmos. Sci.
,
59
,
2647
2664
.
Farrell
,
B. F.
, and
P. J.
Ioannou
,
2002b
:
Optimal perturbation of uncertain systems.
Stochastic Dyn.
,
2
,
395
402
.
Glover
,
K.
,
1984
:
An optimal Hankel-norm approximation of linear multivariable systems and their L-error bounds.
Int. J. Control
,
39
,
1115
1193
.
Högberg
,
M.
,
T. R.
Bewley
, and
D. S.
Henningson
,
2003a
:
Linear feedback control and estimation of transition in plane channel flow.
J. Fluid Mech.
,
481
,
149
175
.
Högberg
,
M.
,
T. R.
Bewley
, and
D. S.
Henningson
,
2003b
:
Relaminarization of Reτ = 100 turbulence using linear state-feedback control.
Phys. Fluids
,
15
,
3572
3575
.
Jovanović
,
M. R.
, and
B.
Bamieh
,
2001
:
Modeling flow statistics using the linearized Navier–Stokes equations. Proc. 40th IEEE Conf. on Decision and Control, Orlando, FL, IEEE, 4944–4949
.
Kalman
,
R. E.
,
1964
:
When is a control system optimal?
J. Basic Eng.
,
83D
,
95
108
.
Kim
,
J.
,
2003
:
Control of turbulence.
Phys. Fluids
,
15
,
1093
1105
.
Kleeman
,
R.
, and
A. M.
Moore
,
1997
:
A theory for the limitation of ENSO predictability due to stochastic atmospheric transients.
J. Atmos. Sci.
,
54
,
753
767
.
Lacarra
,
J.
, and
O.
Talagrand
,
1988
:
Short-range evolution of small perturbations in a barotropic model.
Tellus
,
40A
,
81
95
.
Luenberger
,
D. G.
,
1971
:
An introduction to observers.
IEEE Trans. Autom. Control.
,
16
,
596
602
.
Molteni
,
F.
, and
T. N.
Palmer
,
1993
:
Predictability and finite-time instability of the Northern Hemisphere winter circulation.
Quart. J. Roy. Meteor. Soc.
,
119
,
269
298
.
Palmer
,
T. N.
,
R.
Gelaro
,
J.
Barkmeijer
, and
R.
Buizza
,
1998
:
Singular vectors, metrics, and adaptive observations.
J. Atmos. Sci.
,
55
,
633
653
.
Reddy
,
S.
, and
D. S.
Henningson
,
1993
:
Energy growth in viscous channel flows.
J. Fluid Mech.
,
252
,
209
238
.
Zhou
,
K.
, and
J. C.
Doyle
,
1998
:
Essentials of Robust Control.
Prentice Hall, 411 pp
.

Footnotes

Corresponding author address: Dr. Brian F. Farrell, Harvard University, Pierce Hall 107V, Oxford St. Mail Area H0162, Cambridge, MA 02138. Email: farrell@deas.harvard.edu

1

The reason is that in two-dimensional systems the eigenvalues of the covariance matrix 𝗖 and the stochastic optimal matrix 𝗕 are equal because the 𝗕 and 𝗖 matrices possess two invariants: they have the same trace and also satisfy the relation ℜ [trace (𝗔𝗖)] = ℜ [trace (𝗔𝗕)], which can be easily verified from their respective Lyapunov equations (ℜ denotes the real part). The eigenvalues of 𝗖 and 𝗕 are, in general, different for nonnormal systems with dimension greater than two.