Forward Sensitivity Approach to Dynamic Data Assimilation

Abstract

The least squares fit of observations with known error covariance to a strong-constraint dynamical model has been developed through use of the time evolution of sensitivity functions—the derivatives of model output with respect to the elements of control (initial conditions, boundary conditions, and physical/empirical parameters). Model error is assumed to stem from incorrect specification of the control elements. The optimal corrections to control are found through solution to an inverse problem. Duality between this method and the standard 4D-Var assimilation using adjoint equations has been proved. The paper ends with an illustrative example based on a simplified version of turbulent heat transfer at the sea/air interface.

1. Introduction

Sensitivity function analysis has proved valuable as a mean to both build models and to interpret their output in chemical kinetics (Rabitz et al. [1], Seefeld and Stockwell [2]) and air quality modeling (Russell et al. [3]). Yet, the ubiquitous systematic errors that haunt dynamical prediction cannot be fully understood with sensitivity functions alone. We now include an optimization component that leads to an improved fit of model to observations. The methodology is termed forward sensitivity method (FSM)—a method based on least squares fit of model to data, but where algorithmic structure and correction procedure are linked to the sensitivity functions. In essence, corrections to control (the initial conditions, the boundary conditions, and the physical and empirical parameters) are found through solution to an inverse problem.

In this paper we derive the governing equations for corrections to control and show their equivalence to equations governing the so-called 4D-Var assimilation method (four-dimensional variational method)—least squares fit of model to observations under constraint (LeDimet and Talagrand [4]). Beyond this equivalence, we demonstrate the value of the FSM as a diagnostic tool that can be used to understand the relationship between sensitivity and correction to control.

We begin our investigation by laying down the dynamical framework for the FSM: general form of the governing dynamical model, the type and representation of model error that can identified through the FSM, and the evolution of the sensitivity functions that are central to execution of the FSM. The dual relationship between 4D-Var/adjoint equations is proved. The step-by-step process of assimilating data by FSM is outlined, and we demonstrate its usefulness by application to a simplified air-sea interaction model.

2. Foundation Dynamics for the FSM

We have included a list of mathematical symbols used in this paper. These symbols and associated nomenclature are found in Table 1.

Table 1 

Symbolism and nomenclature.

2.1. Prediction Equations

Let 𝑥(𝑡)∈𝑅^𝑛 denote the state and let 𝛼∈𝑅^𝑝denote the parameters of a deterministic dynamical system, where 𝑥(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡),…,𝑥_𝑛(𝑡))^𝑇 and 𝛼=(𝛼₁,𝛼₂,…,𝛼_𝑝)^𝑇 are column vectors, 𝑛 and 𝑝 are positive integers, 𝑡≥0 denotes the time, and superscript 𝑇 denotes the transpose of the vector or matrix. Let 𝑓∶𝑅^𝑛×𝑅^𝑝×𝑅→𝑅^𝑛 be a mapping, where 𝑓(𝑥,𝛼,𝑡)=(𝑓₁,𝑓₂,…,𝑓_𝑛)^𝑇 with 𝑓_𝑖=𝑓_𝑖(𝑥,𝛼,𝑡) for 1≤𝑖≤𝑛. The vector spaces 𝑅^𝑛 and 𝑅^𝑝 are called the model space and parameter space, respectively.

Consider a dynamical system described by a system of ordinary nonlinear differential equations of the form 𝑑𝑥𝑑𝑡=𝑓(𝑥,𝛼,𝑡), (1a) or in component form 𝑑𝑥_𝑖𝑑𝑡=𝑓_𝑖(𝑥,𝛼,𝑡), (1b)where 𝑑𝑥/𝑑𝑡 denotes the time derivative of the state 𝑥(𝑡), with 𝑥(0)∈𝑅^𝑛 the given initial condition. The control vector for the model is given by 𝑐=(𝑥(0),𝛼)∈𝑅^𝑛×𝑅^𝑝, the combination of initial condition and parameters referred to as the control space. It is tacitly assumed that the map of 𝑓 in (1a) and (1b) is such that the solution 𝑥(𝑡)=𝑥(𝑡,𝑥(0),𝛼)=𝑥(𝑡,𝑐) exists and is unique. It is further assumed that 𝑥(𝑡) has a smooth dependence on the control vector 𝑐 such that the first 𝑘 (≥1) partial derivatives of 𝑥(𝑡) with respect to the components of 𝑐 also exist. The solution 𝑥(𝑡) of (1a) and (1b) is known as the deterministic forecast of the state of the system at time 𝑡>0. If the map 𝑓(⋅) in (1a) and (1b) depends explicitly on 𝑡, then this system is called a time varying or nonautonomous system; if 𝑓(⋅) does not depend on 𝑡, then the system is known as a time invariant or autonomous system.

Let 𝑧(𝑡)∈𝑅^𝑚 be the observation vector obtained from the field measurements at time 𝑡≥0. Let ℎ∶𝑅^𝑛→𝑅^𝑚 be the mapping from the model space 𝑅^𝑛 to the observation space 𝑅^𝑚.

If 𝑥(𝑡) denotes the (unknown) true state, then we assume that 𝑧(𝑡) is given by𝑧(𝑡)=ℎ(𝑥(𝑡))+𝑣(𝑡), (2) where 𝑣(𝑡)∈𝑅^𝑚 is the additive (unobservable and unavoidable) noise. The mapping ℎ(⋅) is known as the forward operator or the observation operator. It is further assumed that 𝑣(𝑡) is a white Gaussian noise with mean zero possessing a known covariance matrix 𝐑(𝑡)∈𝑅^𝑚×𝑚. That is, 𝑣(𝑡)∼𝑁(0,𝐑(𝑡)).

2.2. A Classification of Forecast Errors

The forecast error 𝑒_𝐹(𝑡)∈𝑅^𝑚 is defined as follows:𝑒_𝐹(𝑡)≡𝑧(𝑡)−ℎ(𝑥(𝑡))=𝑏(𝑡)+𝑣(𝑡), (3) the sum of a deterministic part 𝑏(𝑡)=ℎ(𝑥(𝑡))−ℎ(𝑥(𝑡)), (4) and the random part 𝑣(𝑡) induced by the observation noise. Our immediate goal is to analyze and isolate sources and types of forecast errors.

First, if the model map 𝑓(⋅) and the forward operator ℎ(⋅) are without error, that is, exact, and if the control vector 𝑐 is also error free, then the deterministic forecast 𝑥(𝑡) must be correct in the sense that 𝑥(𝑡)=𝑥(𝑡), the true state. Then from (3), the forecast error is purely random or white Gaussian noise. That is, 𝑒_𝐹(𝑡)=𝑣(𝑡). (5)

Second, if 𝑓(⋅) and ℎ(⋅) are known perfectly but 𝑐 has an error, then the forecast 𝑥(𝑡) will have a deterministic error induced by the incorrect control vector. In such a case, we can, in principle, decompose the forecast error as a sum 𝑒_𝐹(𝑡)=𝑏(𝑐,𝑡)+𝑣(𝑡), (6) where the deterministic part, 𝑏(𝑐,𝑡)=ℎ(𝑥(𝑡))−ℎ(𝑥(𝑡)), is purely a function of the control vector error.

Third, if 𝑓(⋅), ℎ(⋅), and 𝑐 are in error, then the forecast error can be represented as 𝑒_𝐹(𝑡)=𝑏(𝑐,𝑓,ℎ,𝑡)+𝑣(𝑡), (7) where the deterministic part 𝑏(𝑐,𝑓,ℎ,𝑡) may have a complex (additive and/or multiplicative) dependence on errors in 𝑐, 𝑓(⋅), and ℎ(⋅).

The following assumption is key to our analysis that follows. The model of choice is faithful to the phenomenon under study. The system is predicted with fidelity—the forecasted state is creditable and useful in understanding the dynamical processes that underpin the phenomenon. Certainly, the forecast will generally exhibit some error, but the primary physical processes are included; that is, the vector field 𝑓includes the pertinent physical processes. In this situation the forecast error stems from erroneously specified elements of control. Thus, in our study the forecast error assumes the form shown in (6). Dee’s work [5] contains a very good discussion of the estimation of the bias 𝑏 in (7) arising from errors in the model and/or observations.

2.3. Dynamics of First-Order Sensitivity Function Evolution

Since our approach is centered on sensitivity functions, we develop the dynamics of evolution of the forward sensitivities in this section.

Differentiating both sides of (1b) with respect to 𝛼_𝑗, interchanging the order of differentiation on the left-hand side, we obtain 𝜕(𝑑𝑥_𝑖/𝑑𝑡)𝜕𝛼_𝑗=𝑑(𝑑𝑡𝜕𝑥_𝑖𝜕𝛼_𝑗)=^𝑛_𝑘=1[(𝜕𝑓_𝑖𝜕𝑥_𝑘)(𝜕𝑥_𝑘𝜕𝛼_𝑗)+𝜕𝑓_𝑖𝜕𝛼_𝑗] (8) for 1≤𝑖≤𝑛 and 1≤𝑗≤𝑝 with 𝜕𝑥_𝑖(0)/𝜕𝛼_𝑗=0 as the initial condition.

These 𝑛𝑝 equations can be succinctly written in matrix form as𝑑𝑑𝑡[𝐷_𝛼(𝑥)]=𝐷_𝑥(𝑓)𝐷_𝛼(𝑥)+𝐷_𝛼(𝑓) (9) with 𝐷_𝛼(𝑥(0))=0 as initial condition. This system of linear time-varying ordinary nonhomogeneous differential equations describe the evolution of the elements of 𝐷_𝛼(𝑥)=[𝜕𝑥_𝑖/𝜕𝛼_𝑗]∈𝑅^𝑛×𝑝, where the Jacobian matrices 𝐷_𝑥(𝑓)∈𝑅^𝑛×𝑛 and 𝐷_𝛼(𝑓)∈𝑅^𝑛×𝑝 are given by𝐷_𝑥(𝑓)=[𝜕𝑓_𝑖𝜕𝑥_𝑗𝐷],_𝛼(𝑓)=[𝜕𝑓_𝑖𝜕𝛼_𝑗]. (10)(11) Similarly, by differentiating both sides of (1b) with respect to 𝑥_𝑗(0), we obtain𝜕(𝑑𝑥_𝑖/𝑑𝑡)𝜕𝑥_𝑗=𝑑(0)(𝑑𝑡𝜕𝑥_𝑖𝜕𝑥_𝑗(0))=^𝑛_𝑘=1(𝜕𝑓_𝑖𝜕𝑥_𝑘)(𝜕𝑥_𝑘𝜕𝑥_𝑗(0)) (12) for 1≤𝑖,𝑗≤𝑛, with 𝜕𝑥_𝑖(0)/𝜕𝑥_𝑗(0)=𝛿_𝑖𝑗, where 𝛿_𝑖𝑗 is the standard Kronecker delta. These 𝑛² equations can be succinctly represented as 𝑑𝑑𝑡[𝐷_𝑥(0)(𝑥)]=𝐷_𝑥(𝑓)𝐷_𝑥(0)(𝑥) (13) with 𝐷_𝑥(0)(𝑥(0))=𝐼, the identity matrix. This system of linear, time-varying homogeneous equations governs the evolution of the elements of the matrix 𝐷_𝑥(0)(𝑥)=[𝜕𝑥_𝑖/𝜕𝑥_𝑗(0)]∈𝑅^𝑛×𝑝. Notice that (9) and (13) are independent of the observations and have the same system matrix 𝐷_𝑥(𝑓) on the right-hand sides; thus, the homogeneous solutions to (9) and (13) have the same structure.

The evolution of the sensitivities (solution to (9) and (13)) is dependent on the solution to the governing dynamical equations ((1a) and (1b)). Generally, these equations are solved numerically using the standard fourth-order Runge-Kutta method. Rabitz et al. work [1] contains more details relating to solutions of (9) and (13). In special cases such as in air quality modeling, the sensitivity equations (9) and (13) exhibit extreme stiffness. Special methods are needed to handle the inherent stiffness of these equations. Seefeld and Stockwell work [2] includes a discussion of these issues. Gear’s work [6] is a good reference for a general discussion of stiff equations.

3. Duality between the FSM and 4D-Var Based on Adjoint Method

Let 𝑧(𝑡₁),𝑧(𝑡₂),…,𝑧(𝑡_𝑁) be a sequence of 𝑁 observation vectors at times 𝑡₀=0<𝑡₁<𝑡₂⋯<𝑡_𝑁. The goal is to use these observations to improve the estimate of the control vector 𝑐. This estimation problem is recast as a constrained minimization of an objective function 𝐽∶𝑅^𝑛×𝑅^𝑝→𝑅 given by1𝐽(𝑐)=2^𝑁𝑖=1[𝑧(𝑡_𝑖)−ℎ(𝑥(𝑡_𝑖))]^𝑇𝐑⁻¹(𝑡_𝑖)[𝑧(𝑡_𝑖)−ℎ(𝑥(𝑡_𝑖))], (14) where the model state 𝑥(𝑡) evolves according to (1a), (1b), and 𝐑(𝑡_𝑖)∈𝑅^𝑚×𝑚 is the known covariance of the observational errors 𝑣(𝑡_𝑖) at time 𝑡_𝑖,1≤𝑖≤𝑁.

Fundamental to minimizing (14) is the computation of the gradient of 𝐽(𝑐), denoted by ∇_𝑐𝐽(𝑐). In the following we describe two ways of characterizing ∇_𝑐∇𝐽(𝑐)=⌊^𝑥(0)𝐽∇−−−_𝛼𝐽⌋, (15) where ∇_𝑥(0)𝐽∈𝑅^𝑛 and ∇_𝛼𝐽∈𝑅^𝑝.

3.1. The Adjoint Method

This method is based on the basic principle that if 𝛿𝐽 is the first variation of 𝐽(𝑐) induced by the perturbation 𝛿𝑐 in 𝑐, then 𝛿𝐽=⟨∇_𝑐𝐽,𝛿𝑐⟩, (16) where ⟨𝑎,𝑏⟩ denotes the standard inner product of two vectors 𝑎 and 𝑏 of the same dimension. Once the first variation 𝛿𝐽 is determined, the gradient can be found. In the following we exploit two basic properties of inner product:

(i) linearity ⟨𝑎+𝑏,𝑥⟩=⟨𝑎,𝑥⟩+⟨𝑏,𝑥⟩, (17)

(ii) adjoint property ⟨𝑎,𝐺𝑏⟩=⟨𝐺^𝑇𝑎,𝑏⟩, (18) where 𝐺^𝑇 is the transpose or the adjoint of the matrix 𝐺.

From first principles (Chapter 24 in Lewis et al. [7]), it can be verified that the first variation 𝛿𝐽 of 𝐽 in (14) is given by 𝛿𝐽=−^𝑁_𝑘=1⟨𝐑⁻¹(𝑡_𝑘)𝑒_𝐹(𝑡_𝑘),𝐷_𝑥(ℎ)𝛿𝑥(𝑡_𝑘)⟩, (19) where the forecast error given by 𝑒_𝐹(𝑡_𝑘)=𝑧(𝑡_𝑘)−ℎ(𝑥(𝑡_𝑘𝐷)),_𝑥(ℎ)=[𝜕ℎ_𝑖𝜕𝑥_𝑗]∈𝑅^𝑚×𝑛 (20) is the Jacobian of forward operator ℎ(𝑥) with respect to 𝑥.

By invoking the adjoint property in (17), (19) becomes 𝛿𝐽=−^𝑁_𝑘=1⟨𝜂(𝑡_𝑘),𝛿𝑥(𝑡_𝑘)⟩, (21) where 𝜂(𝑡_𝑘)=𝐷^𝑇_𝑥(ℎ)𝐑⁻¹(𝑡_𝑘)𝑒_𝐹(𝑡_𝑘). (22) The first variation 𝛿𝑥(𝑡_𝑘) in 𝑥(𝑡) at 𝑡=𝑡_𝑘 resulting from the perturbation 𝛿𝑐 in 𝑐 is given by (A.4) in Appendix A. Using (A.4) and the linearity of the inner product it follows that 𝛿𝐽=−^𝑁_𝑘=1⟨𝜂(𝑡_𝑘),𝑀(𝑡_𝑘−,0)𝛿𝑥(0)⟩_𝑁𝑘=1⟨𝜂(𝑡_𝑘),𝑀(𝑡_𝑘,0)𝐿(𝑡_𝑘,0)𝛿𝛼⟩, (23) where we will refer to the first term on the right-hand side of (23) as “Term I” and the second as “Term II.” Using the adjoint and linearity property in (17), we getTermI=−⟨^𝑁_𝑘=1𝑀^𝑇(𝑡_𝑘,0)𝜂(𝑡_𝑘),𝛿𝑥(0)⟩ (24) from which we obtain∇_𝑥(0)𝐽=−^𝑁_𝑘=1𝑀^𝑇(𝑡_𝑘,0)𝜂(𝑡_𝑘). (25) Similarly, rewriting Term II as TermII=−⟨^𝑁_𝑘=1𝐿^𝑇(𝑡_𝑘,0)𝑀^𝑇(𝑡_𝑘,0)𝜂(𝑡_𝑘),𝛿𝛼⟩, (26) we get ∇_𝛼𝐽=−^𝑁_𝑘=1𝐿^𝑇(𝑡_𝑘,0)𝑀^𝑇(𝑡_𝑘,0)𝜂(𝑡_𝑘). (27) Hence, we obtain the components of ∇_𝑐𝐽 which are used in conjunction with the minimization algorithm to find the optimal 𝑐^∗ that minimizes 𝐽(𝑐) in (14).

We conclude this discussion with an efficient recursive method for evaluating the expressions in (25) and (27). Define𝜆_𝑁=𝜂(𝑡_𝑁), (28) and for 𝑘=𝑁−1,𝑁−2,…,2,1,0,𝜆_𝑘=𝑀^𝑇(𝑡_𝑘+1,𝑡_𝑘)𝜆_𝑘+1+𝜂(𝑡_𝑘). (29) It can be verified that ∇_𝑥(0)𝐽=𝑀^𝑇(𝑡₁,0)𝜆₁ in (25).

Details on the recursive computation of (27) are given in Appendix C.

3.2. Sensitivity-Based Approach

Let us first consider the special case when 𝑁=1. Then (14) becomes𝐽₁(𝑐⁻1)=2[𝑧−ℎ(𝑥)]^𝑇𝐑⁻¹(𝑡)[𝑧−ℎ(𝑥)], (30) where we recall that 𝑐∈𝑅^𝑛×𝑅^𝑝 and 𝑥=𝑥(𝑡,𝑐) is the solution of the model equations (1a) and (1b) at time 𝑡.

Setting 𝑞=(𝑛+𝑝) and 𝐴=𝐑⁻¹(𝑡), the expression 𝐽₁(𝑐) in (30) becomes identical to 𝑄(𝑐) in (B.16) (Appendix B). Hence, by using (B.20) it follows that ∇_𝑐𝐽₁(𝑐)=−𝐷^𝑇_𝑐(𝑥)𝜂(𝑡), (31) where 𝜂(𝑡) is given by (22) and 𝐷^𝑇_𝑐𝐷(𝑥)=⌊^𝑇𝑥(0)𝐷(𝑥)−−−−−_𝑇𝛼(𝑥)⌋∈𝑅^{(𝑛+𝑝)×𝑛}. (32) Now comparing (32) with (25)–(27), we obtain the duality relation:𝐷_𝑥(0)𝐷(𝑥)⟷𝑀(𝑡,0),_𝛼(𝑥)⟷𝑀(𝑡,0)𝐿(𝑡,0). (33) The sensitivity matrices on the left-hand side of (33) are obtained by solving the forward sensitivity equations, whereas the matrix products 𝑀(𝑡,0) and 𝐿(𝑡,0) are obtained by integrating along the path as given by (A.3) and (A.5) (Appendix A).

The primary advantage of the sensitivity-based approach is that it provides a natural interpretation of the expression for the gradient in (31). Recall from (22) that 𝜂(𝑡)∈𝑅^𝑛 is the weighted version of the forecast error 𝐑⁻¹(𝑡)𝑒_𝐹(𝑡)∈𝑅^𝑚 mapped onto the model space by the linear operator 𝐷^𝑇_𝑥(ℎ). Also the 𝑖th column of 𝐷^𝑇_𝑐(𝑥) is the sensitivity of the 𝑖th component 𝑥_𝑖(𝑡) of the state vector 𝑥(𝑡)∈𝑅^𝑛 with respect to the control vector 𝑐∈𝑅^𝑞 given by (𝜕𝑥_𝑖/𝜕𝑐₁,𝜕𝑥_𝑖/𝜕𝑐₂,…,𝜕𝑥_𝑖/𝜕𝑐_𝑞). Thus, it follows from (31) that ∇_𝑐𝐽₁ is a linear combination of the columns of 𝐷^𝑇_𝑐(𝑥) which are the sensitivity vectors, where the coefficients of the linear terms are the components of 𝜂(𝑡). Thus, columns of 𝐷^𝑇_𝑐(𝑥) with large norms that are associated with large forecast errors will be dominant in determination of the components of ∇𝐽₁(𝑐). In other words, we gain some insight into the interplay between corrections to control and forecast errors—something that can be seen through a careful examination of the sensitivity vector at various times from initial state to forecast horizon. (The illustrative example in Section 5 further explores this diagnostic function.) Expression (31) also enables us to isolate the effect of different components 𝑥_𝑖 of 𝑥(𝑡) on the performance index 𝐽₁(𝑐).

For the general case of observations at multiple times, (31) assumes the following form:∇_𝑐𝐽=−^𝑁_𝑖=1𝐷^𝑇_𝑐(𝑥(𝑡_𝑖))𝜂(𝑡_𝑖). (34) This gradient is the sum of linear combinations of the columns of 𝐷^𝑇_𝑐(𝑥(𝑡_𝑖)) at various time instances. With so many directions (the directions associated with the columns of 𝐷^𝑇_𝑐(𝑥(𝑡_𝑖)) contributing to the components of ∇_𝑐𝐽, the connection between sensitivity, and forecast error is obscured.

4. Data Assimilation Using Sensitivity

We seek to find the solution to the following problem using the FSM. Given 𝑓(⋅) and ℎ(⋅), the control vector 𝑐, the observation 𝑧(𝑡), and the error covariance of observations 𝐑(t), find a correction 𝛿𝑐 to the control vector such that the new model forecast starting from (𝑐+𝛿𝑐) will render the forecast error 𝑒_𝐹(𝑡) purely random; that is, the systematic forecast error is removed and accordingly 𝐸(𝑒_𝐹(𝑡))=0.

We start by quantifying the change Δ𝑥 in the solution 𝑥(𝑡)=𝑥(𝑡,𝑐) resulting from a change 𝛿𝑐 in 𝑐. Invoking the standard Taylor series expansion, we obtainΔ𝑥=𝑥(𝑡,𝑐+𝛿𝑐)−𝑥(𝑡,𝑐)=_𝑗𝛿^𝑗𝑥, (35) where 𝛿^𝑘𝑥 is the 𝑘th variation of 𝑥(𝑡), the fraction of the total change that can be accounted by the 𝑘th partial derivatives of 𝑥(𝑡) with respect to 𝑐 and the perturbation 𝛿𝑐. Since practical considerations dictate that the total number of correction terms on the right-hand side of (35) be finite, we often settle for an approximation of only 𝑘 terms (𝑘 generally ≤ 2). This is a tradeoff between improved accuracy resulting from a large value of 𝑘 and the complexity of the resulting inverse problem. Although we have developed the methodology for second-order analysis (𝑘 = 2, where Δ𝑥 is approximated by the sum 𝛿𝑥+𝛿²𝑥) (Lakshmivarahan and Lewis [8]), our development will follow the first-order analysis (𝑘=1, where Δ𝑥 is given by the first variation 𝛿𝑥). Second-order analysis is justified when 𝛿²𝑥 is a significant fraction of 𝛿𝑥—this occurs when 𝑓(𝑥) and/or ℎ(𝑥) exhibit strong nonlinearity. It is further shown in Section 5 that iterative application of the first-order method often leads to improved results.

4.1. First-Order Analysis

From first principles and (35) we obtain Δ𝑥≈𝛿𝑥=𝐷_𝑥(0)(𝑥)𝛿𝑥(0)+𝐷_𝛼(𝑥)𝛿𝛼, (36) where 𝐷_𝑥(0)(𝑥)∈𝑅^𝑛×𝑛 is the Jacobian of 𝑥(𝑡) with respect to 𝑥(0), and 𝐷_𝛼(𝑥)∈𝑅^𝑛×𝑝 is the Jacobian of 𝑥(𝑡) with respect to 𝛼.The matrices 𝐷_𝑥(0)(𝑥) and 𝐷_𝛼(𝑥), found as solutions of (13) and (9), respectively, are known as the first-order sensitivity of the solution 𝑥(𝑡) with respect to 𝑥(0) and 𝛼, respectively, and the elements of these matrices are called sensitivity functions.

4.1.1. Observations at One Time Only

We first consider the case where observation 𝑧(𝑡)∈𝑅^𝑚 is available at one time 𝑡. The first variation 𝛿𝑥 in 𝑥(𝑡) induces a variation Δℎ in the forward operator ℎ(𝑥(𝑡)). Again, by approximating Δℎ by the first variation, we getΔℎ≈𝛿ℎ=𝐷_𝑥(ℎ)𝛿𝑥, (37) where 𝐷_𝑥(ℎ)∈𝑅^𝑚×𝑛 is the Jacobian of ℎ(⋅) with respect to 𝑥 and is given by𝐷_𝑥(ℎ)=[𝜕ℎ_𝑖𝜕𝑥_𝑗],1≤𝑖≤𝑚,1≤𝑗≤𝑛. (38) substituting (36) into (37), we obtain𝑑ℎ=𝐻₁𝑑𝑥(0)+𝐻₂𝑑𝛼, (39) where 𝐻₁(𝑡)=𝐷_𝑥(ℎ)𝐷_𝑥(0)(𝑥)∈𝑅^𝑚×𝑛 and 𝐻₂(𝑡)=𝐷_𝑥(ℎ)𝐷_𝛼(𝑥)∈𝑅^𝑚×𝑝. Setting 𝐻(𝑡)=[𝐻₁(𝑡),𝐻₂(𝑡)]∈𝑅^{𝑚×(𝑛+𝑝)} and 𝜍=(𝜍₁,𝜍₂)∈𝑅^𝑛+𝑝, where 𝜍₁=𝛿𝑥(0) and 𝜍₂=𝛿𝛼, (39) becomes𝑑ℎ=𝐻(𝑡)𝜍. (40) Given the operating point 𝑐, our goal is to find the perturbation 𝛿𝑐 such that the observation is equal to the model counterpart, that is, 𝑧(𝑡)=ℎ(𝑥(𝑡,𝑐+𝛿𝑐))≈ℎ(𝑥(𝑡))+𝛿ℎ, (41) or𝑒_𝐹(𝑡)≡𝑧(𝑡)−ℎ(𝑥(𝑡))=𝛿ℎ. (42) From (43), it follows that the required perturbation 𝜍∈𝑅^𝑛+𝑝 is obtained by solving the linear inverse problem𝐻(𝑡)𝜍=𝑒_𝐹(𝑡), (43) where 𝐻(𝑡)∈𝑅^{𝑚×(𝑛+𝑝)} and 𝑒_𝐹(𝑡)∈𝑅^𝑚.

4.1.2. Observations at Multiple Times

The above analysis can be readily extended to the case where observations are available at 𝑁 times. We denote these sets of observation vectors by 𝑧(𝑡₁),𝑧(𝑡₂),…,𝑧(𝑡_𝑁), where 0<𝑡₁<𝑡₂<⋯<𝑡_𝑁. The forecast error 𝑒_𝐹(𝑡_𝑖) is given by 𝑒_𝐹(𝑡_𝑖)=𝑧(𝑡_𝑖)−ℎ(𝑥(𝑡_𝑖)),1≤𝑖≤𝑁. (44) Define𝐻₁(𝑡_𝑖)=𝐷_{𝑥(𝑡𝑖)}(ℎ)𝐷_𝑥(0)(𝑥(𝑡_𝑖))∈𝑅^𝑚×𝑛,𝐻₂(𝑡_𝑖)=𝐷_{𝑥(𝑡𝑖)}(ℎ)𝐷_𝛼(𝑥(𝑡_𝑖))∈𝑅^𝑚×𝑝. (45) Then at time 𝑡_𝑖 we have𝐻(𝑡_𝑖)𝜍=𝑒_𝐹(𝑡_𝑖),1≤𝑖≤𝑁, (46) where 𝐻(𝑡_𝑖)=[𝐻₁(𝑡_𝑖),𝐻₂(𝑡_𝑖)]∈𝑅^{𝑚×(𝑛+𝑝)}. (47) Now define a matrix 𝐻∈𝑅^{𝑁𝑚×(𝑛+𝑝)} and a vector 𝑒_𝐹∈𝑅^𝑁𝑚 as𝐻=𝐻(𝑡₁)𝐻(𝑡₂)⋮𝐻(𝑡_𝑁),𝑒_𝐹=𝑒_𝐹(𝑡₁)𝑒_𝐹(𝑡₂)⋮𝑒_𝐹(𝑡_𝑁). (48) Then, the N relations in (46) can be succinctly denoted by𝐻𝜍=𝑒_𝐹. (49) A number of special cases arise depending on (a) the value of 𝑁𝑚 relative to (𝑛+𝑝), namely, over (under) determined cases when 𝑁𝑚>(<)(𝑛+𝑝) and (b) the rank of the matrix 𝐻(𝑡), namely, 𝐻(𝑡) is of full rank or rank deficient. In all these cases, the linear inverse problem (43) is recast as a minimization problem using the standard least squares framework (Lawson and Hanson [9]). The resulting minimization problem can then be solved using one of many standard methods, for example, the conjugate gradient method (Lewis et al. [7]; Nash and Sofer [10]).

As an illustration, consider the case when 𝑁𝑚>(𝑛+𝑝) and that the rank of 𝐻 is (𝑛+𝑝), that is, full rank. The solution 𝜍 is then obtained by minimizing the weighted least squares criterion𝐽_𝑁1(𝜍)=2^𝑁𝑖=1(𝑒_𝐹(𝑡_𝑖)−𝐻(𝑡_𝑖)𝜍)^𝑇𝐑⁻¹(𝑡_𝑖)(𝑒_𝐹(𝑡_𝑖)−𝐻(𝑡_𝑖=1)𝜍)2(𝑒_𝐹−𝐻𝜍)^𝑇𝐑⁻¹(𝑒_𝐹−𝐻𝜍), (50) where 𝐑=Diag[𝐑(𝑡_i),𝐑(𝑡₂),…,𝐑(𝑡_N)] (51) is an 𝑁𝑚×𝑁𝑚 diagonal matrix with 𝐑(t_i) as its 𝑖th diagonal block.

Although it is computationally efficient to minimize (50) by using a method like conjugate gradient, there is an advantage to analyze the properties of the optimal solution via the classical approach, that is, by setting the gradient of 𝐽_𝑁(𝜍) to zero. It can be verified that the minimizing 𝐽_𝑁 is found by solving the symmetric linear system[^𝑁_𝑖=1𝐻^𝑇(𝑡_𝑖)𝐑⁻¹(𝑡_𝑖)𝐻(𝑡_𝑖)]𝜍=^𝑁_𝑖=1(𝐻^𝑇(𝑡_𝑖)𝐑⁻¹(𝑡_𝑖)𝑒_𝐹(𝑡_𝑖)) (52) or succinctly as𝜍_LS=(𝐻^𝑇𝐑⁻¹𝐻)⁻¹𝐻^𝑇𝐑⁻¹𝑒_𝐹, (53) where 𝐻,𝑒_𝐹, and 𝐑 are defined in (48) and (51), and subscript “LS” refers to the least squares solution.

From the discussion relating to the classification of forecast errors, recall that the forecast error inherits its randomness from the (unobservable) observation noise. The vector 𝑒_𝐹 on the right hand side of (53) is random and hence the solution 𝜍 of (53) is also random.

Since we are interested in forecast errors in response to incorrect control, we have 𝑒_𝐹(𝑡_𝑖)=𝑏(𝑐,𝑡_𝑖)+𝑣(𝑡_𝑖)for1≤𝑖≤𝑁. (54) Now define𝑏=𝑏(𝑐,𝑡₁)𝑏(𝑐,𝑡₂)⋮𝑏(𝑐,𝑡_𝑁),𝑣=𝑣(𝑡₁)𝑣(𝑡₂)⋮𝑣(𝑡_𝑁). (55) Hence, the vector 𝑒_𝐹 in (48) can be expressed as𝑒_𝐹=𝑏+𝑣 (56) with 𝐸(𝑒_𝐹)=𝑏 since 𝐸(𝑣)=0. Substituting (56) into (53) and taking the expectation give𝐸(𝜍_LS)=[(𝐻^𝑇𝐑⁻¹𝐻)⁻¹𝐻^𝑇𝐑⁻¹]𝑏. (57) Thus, the expected value of the correction to control is indeed a linear function of the forecast error 𝑏 itself. It can be verified (Lewis et al. [7]) that the covariance of the least squares estimate [(53)] is given by cov(𝜍_LS)=(𝐻^𝑇𝐑⁻¹𝐻)⁻¹=[∇²𝐽_𝑁(𝜍)]⁻¹, (58) where ∇²𝐽_𝑁(𝜍) is the Hessian of 𝐽_𝑁(𝜍) in (50).

5. A Practical Example: Air/Sea Interaction

We choose a simple but nontrivial differential equation to demonstrate the applicability of the forward-sensitivity method to identification of error in a dynamical model. We break this discussion into three parts as follows: (1) the model, (2) discussion of the diagnostic value of FSM, and (3) numerical experiments with data assimilation using FSM.

5.1. The Model

Consider the case where cold continental air moves out over an ocean of constant surface temperature. We follow a column of air in a Lagrangian frame; that is, the column of air moves with the prevailing low-level wind speed. Turbulent transfer of heat from the ocean to air warms the column. The governing equation is taken to be𝑑𝜃=𝐶𝑑𝜏^𝑇𝑉𝐻(𝜃_𝑠−𝜃), (59) where 𝜃: temperature of the air column (^∘C),𝜃_𝑠: temperature of the sea surface (^∘C),𝐶_𝑇: turbulent heat exchange coefficient (nondimensional),𝑉: speed of air column (ms⁻¹),𝐻: height of the column (mixed layer)(m),𝜏: time (h).Equation (59) is nondimensionalized by the following scaling:𝜃=𝜋⋅𝑥(𝜋=1.0^∘𝜃C),_𝑠=𝜋⋅𝑥_𝑠,𝜏=𝑇⋅𝑡(𝑇=1.0h). (60) The governing equation then takes the form𝑑𝑥(𝑡)=𝐶𝑑𝑡^𝑇𝑉𝑇𝐻(𝑥_𝑠−𝑥(𝑡)). (61) Assuming 𝐻∼150m, 𝑉∼10ms⁻¹, 𝐶_𝑇∼10⁻³, then𝐶𝑘=^𝑇𝑉𝑇𝐻≈14. (62) Thus, we take our governing equation to be𝑑𝑥(𝑡)𝑑𝑡=𝑓(𝑥,𝛼)=𝑘(𝑥_𝑠−𝑥(𝑡)), (63) where 𝑘=0.25. The solution to (63) is𝑥(𝑡,𝑥(0),𝛼)=[𝑥(0)−𝑥_𝑠]𝑒^−𝑘𝑡+𝑥_𝑠 (64) with 𝑐=(𝑥(0),𝛼)∈𝑅×𝑅², and 𝛼=(𝑥_𝑠,𝑘)^𝑇∈𝑅², that is, 𝑛=1 and 𝑝 = 2. The solution depends linearly on 𝑥(0) and 𝑥_𝑠 but nonlinearly on 𝑘.

There are three elements of control: initial condition, 𝑥(0), boundary condition, 𝑥_𝑠, and parameter, 𝑘.

5.2. Diagnostic Aspects of FSM

The Jacobians of 𝑓 with respect to 𝑥 and 𝛼 are given by 𝐷_𝑥𝐷(𝑓)=−𝑘∈𝑅,_𝛼(𝑓)=[𝑘,𝑥_𝑠−𝑥(𝑡)]∈𝑅^1×2, (65) and the Jacobians of the solution 𝑥(𝑡) with respect to 𝛼 and 𝑥(0) are given by𝐷_𝛼(𝑥(𝑡))=[𝜕𝑥(𝑡)𝜕𝑥_𝑠,𝜕𝑥(𝑡)𝜕𝑘]∈𝑅^1×2,𝐷_𝑥(0)(𝑥(𝑡))=[𝜕𝑥(𝑡)𝜕𝑥(0)]∈𝑅. (66) From (9) and (13) the evolution of the forward sensitivities is given by 𝑑(𝑑𝑡𝜕𝑥(𝑡)𝜕𝑥_𝑠)=−𝑘(𝜕𝑥(𝑡)𝜕𝑥_𝑠𝑑)+𝑘,(𝑑𝑡𝜕𝑥(𝑡)𝜕𝑘)=−𝑘(𝜕𝑥(𝑡)𝜕𝑘)+(𝑥_𝑠𝑑−𝑥(𝑡)),(𝑑𝑡𝜕𝑥(𝑡)𝜕𝑥(0))=−𝑘(𝜕𝑥(𝑡)𝜕𝑥(0)), (67)(68)(69) where [𝜕𝑥(𝑡)𝜕𝑥_𝑠]_𝑡=0=0=[𝜕𝑥(𝑡)]𝜕𝑘_𝑡=0,[𝜕𝑥(𝑡)]𝜕𝑥(0)_𝑡=0=1. (70) Either by solving (67)–(70) or by computing directly from (64), it can be verified that the required sensitivities evolve according to𝜕𝑥(𝑡)𝜕𝑥(0)=𝑒^−𝑘𝑡,𝜕𝑥(𝑡)𝜕𝑥_𝑠=1−𝑒^−𝑘𝑡,𝜕𝑥(𝑡)𝜕𝑘=[𝑥_𝑠−𝑥(0)]𝑡𝑒^−𝑘𝑡. (71) The plots of the solution and the three sensitivities are given in Figure 1.

(a) Solution 𝑥(𝑡)

(b) Sensitivity of 𝑥(𝑡) w.r.t. 𝑥(0)

(c) Sensitivity of 𝑥(𝑡) w.r.t. 𝑘

(d) Sensitivity of 𝑥(𝑡) w. r. t 𝑥𝑠

Figure 1 

Evolution of the solution of 𝑥(𝑡) and its sensitivities to 𝑐=(𝑥(0),𝑥𝑠,𝑘)=(1.0,11.0,0.25).

Let 𝑧(𝑡) be the direct observation of the state 𝑥(𝑡), namely, 𝑧(𝑡)=𝑥(𝑡)+𝑣(𝑡). (72) In this case, ℎ(𝑡) is the forecast variable and therefore 𝐷_𝑥(ℎ)=1. Then𝑒_𝐹(𝑡)=𝑧(𝑡)−𝑥(𝑡) (73) is the forecast error.

Following the developments in Section 4 for the case of a single observation described by (39)–(46), we obtain the analog of (46) as𝐻𝜍=𝑒_𝐹(𝑡), (74) where 𝐻=𝐻(𝑡)=[𝜕𝑥(𝑡)/𝜕𝑥(0),𝜕𝑥(𝑡)/𝜕𝑥_𝑠,𝜕𝑥(𝑡)/𝜕𝑘] is the forward sensitivity vector and 𝜍=[𝛿𝑥(0),𝛿𝑥_𝑠,𝛿𝑘]^𝑇. Clearly, (74) corresponds to an under-determined linear least squares problem, whose optimal solution is given by [7, chapter 5]𝜍_LS=𝐻^𝑇(𝐻𝐻^𝑇)⁻¹𝑒_𝐹𝐻(𝑡)=^𝑇𝑒‖𝐻‖^𝐹(𝑡)‖𝐻‖, (75) where 𝐻^𝑇/‖𝐻‖ is the unit forward sensitivity vector at time 𝑡 and 𝑒_𝐹(𝑡)/‖𝐻‖ is a scalar that is the forecast error normalized by the norm of the forward sensitivity vector 𝐻. In other words, the direction of the optimal corrections to the control that annihilate the forecast error is a constant multiple of the unit forward sensitivity vector.

If we assume the following control vector: 𝑐=(1.0,11.0,0.25), that is, [𝑥(0)=1^∘C, 𝑥_𝑠=11^∘C, and 𝑘=0.25 (nondimensional)], we get 𝐻^𝑇=𝑒^−0.25𝑡1−𝑒^−0.25𝑡10𝑡𝑒^−0.25𝑡=𝜕𝑥(𝑡)𝜕𝑥(0)𝜕𝑥(𝑡)𝜕𝑥_𝑠𝜕𝑥(𝑡)𝜕𝑘. (76) The time variations of elements of 𝐻 are given in Table 2 (also refer to Figure 1). From this table, it is clear that the direction of corrections to control varies as 𝑡 increases. At 𝑡=0, the corrections lie in the direction (1,0,0)^𝑇, where 𝑥(𝑡) is only sensitive to the initial condition 𝑥(0). For large 𝑡, the corrections lie in the direction (0,1,0)^𝑇, where 𝑥(𝑡) is only sensitive to the boundary condition 𝑥_𝑠 (the sea-surface temperature). For intermediate times, all the components of control have nonzero sensitivities. 𝜕𝑥(𝑡)/𝜕𝑘 reaches its maximum at 𝑡=4.0.

Table 2 

Sensitivities of model state to initial conditions and parameters for the sea/air turbulent transfer model.

5.3. Numerical Experiments

We assume that the incorrect control vector is 𝑐^′=(𝑥^′(0),𝑥^′_𝑠,𝑘^′)=(2.0,10.0,0.30); in dimensional form, 𝑥^′(0)=2^∘C, 𝑥^′_𝑠=10^∘C, and 𝑘^′=0.30 (non-dimensional). Thus, for an ideal correction to control,𝛿𝑐=(−1.0,+1.0,−0.05). (77) To mimic reality, the correction process uses sensitivity functions that stem from the erroneous solution, that is, where the incorrect control is used.

We have explored both the goodness and failure of recovery of control under two different scenarios, where either 3 or 6 observations are used to recover the control vector. Since there are 3 unknowns, the case for 6 observations is an over-determined system.

We execute numerical experiments where the observations are spread over different segments of the temporal range—generally divided into an “early stage” and a “saturated stage.” By saturated stage we refer to that segment where the solution becomes close to the asymptotic state, that is, 𝑥→𝑥_𝑠. The dividing time between these segments is arbitrary; but generally, based on experimental results to follow, we divide the segments at 𝑡 = 24 where 𝑥(24) = 10.975, 0.025 less than 𝑥_𝑠 = 11.0 (see Figure 1).

The following general statement can be made. If more than one of the observations falls in the saturated zone, the matrix becomes illconditioned. As can be seen from (39) and the plots of sensitivity functions in Figure 1, 𝜕𝑥/𝜕𝑥_𝑠→1 and 𝜕𝑥/𝜕𝑥(0) and 𝜕𝑥/𝜕𝑘→ 0 as 𝑡→∞. Accordingly, if two of the observations are made in the saturated zone, this induces linear dependency between the associated rows of the 𝐻 matrix and in turn leads to ill-conditioning. This illconditioning is exhibited by a large value of the condition number, the ratio of the largest to the smallest eigenvalue of the matrix 𝐻^𝑇𝐻. The inversion of this matrix is central to the optimal adjustments of control (see (55)).

Illconditioning can also occur as a function of the observation spacing in certain temporal segments. This is linked to lack of variability or lack of change of sensitivity from one observation time to another. And, as can be seen in Figure 1, the absolute value of the slope in sensitivity function curves is generally large at the early stages of evolution and small at later stages. As an example, we find satisfactory recovery, 𝛿𝑐= (–0.882, – 0.067, +0.922), when the observations are located at 5.0, 5.1, and 5.2 (a uniform spacing of Δ𝑡= 0.1). Yet, near the saturated state, at 𝑡 = 20.0, 20.1, and 20.2, again a spacing of 0.1, the recovery is poor with the result 𝛿𝑐 = (+5.317, – 0.142, +0.998). The associated condition numbers for these two experiments are 1.0×10³ and 1.0×10⁶, respectively. Similar results follow from the case where 6 observations are taken. In all of these cases, the key factor is the condition number of 𝐻^𝑇𝐻. For our dynamical constraint, a condition number less than ~10⁴ portends a good result.

For the case where we have 6 observations at 𝑡 = 2, 7, 12, 17, 22, and 27, with a random error of 0.01 (standard deviation), we have executed an ensemble experiment with 100 members to recover control. In this case, the condition number is 2.4×10³. Results are plotted three dimensionally and in two-dimensional planes in the space of control, that is, plots of correction in the 𝑥_𝑠/𝑥(0), 𝑥_𝑠/𝑘, and 𝑥(0)/𝑘 planes. Results are shown in Figure 2.

(a) 3D clusterplot

(b) Cluster of corrections of 𝑥(0) versus 𝑘

(c) Cluster of corrections of 𝑥(0) versus 𝑥𝑠

(d) Cluster of corrections of 𝑘 versus 𝑥𝑠

Figure 2 

3D cluster and its projections where corrections are computed using observations at the following times: 𝑡 = 2, 7, 12, 17, 22, and 27. An ensemble of 100 members is used where the standard deviation of the observation noise is 0.01.

Finally, we explore the iterative process of finding corrections to control. Here, the results from the 1st iteration are used to find the new control vector. This vector is then used to make another forecast and find a new set of sensitivity functions. The error of the forecast is obtained, and along with the new sensitivity functions, a second correction to control is found, and so forth. For the experiment with 6 observations that has been discussed in the previous paragraph, we apply this iterative methodology. As can be seen in Figure 3, the correct value of control is found in 3 iterations.

(a) Iterative correction

(b) Projection onto 𝑥(0)-𝑘 plane

(c) Projection onto 𝑘-𝑥𝑠 plane

(d) Projection on to 𝑥(0)-𝑥𝑠 plane

Figure 3 

An illustration of the progression of first-order iterative corrections using six observations at 𝑡 = 2, 7, 12, 17, 22, and 27.

6. Concluding Remarks

The basic contributions of this paper may be stated as follows.

(1) While the 4D-Var has been the primary methodology for operational data assimilation in meteorology/oceanography (Lewis and Lakshmivarahan [11]), and while forward sensitivity has been a primary tool for reaction kinetics and chemistry (Rabitz et al. [1]) and air quality modeling (Russell et al. [3], to our knowledge these two methodologies have not been linked. We have now shown that the method of computing the gradient of 𝐽(𝑐) by these two approaches exhibits a duality hitherto unknown.

(2) By treating the forward sensitivity problem as an inverse problem in data assimilation, we are able to understand the fine structure of the forecast error. This is not possible with the standard 4D-Var formulation using adjoint equations.

(3) While it is true that computation of the evolution of the forward sensitivity involves computational demands beyond those required for solving the adjoint equations in the standard 4D-Var methodology, there is a richness or augmentation to the information that comes with this added computational demand. In essence, it allows us to make judicious decisions on placement of observations through understanding of the time dependence of correction to control.

Appendices

A. Dynamics of Evolution of Perturbations

Let 𝛿𝑐=(𝛿𝑥^𝑇(0),𝛿𝛼^𝑇)^𝑇∈𝑅^𝑛×𝑅^𝑝 be the perturbation in the control vector 𝑐 and 𝛿𝑥(𝑡) the resulting perturbation in the state 𝑥(𝑡) induced by the dynamics (1a) and (1b). Our goal is to derive the dynamics of evolution of 𝛿𝑥(𝑡).

From first principles, the evolutionary dynamics of 𝛿𝑥(𝑡) are given by the variational equation (Hirsch et al. [12]) or the tangent linear model (Lewis et al. [7])𝑑𝑑𝑡𝛿𝑥(𝑡)=𝐴(𝑡)𝛿𝑥(𝑡)+𝐵(𝑡)𝛿𝛼, (A.1) where the Jacobians 𝐴(𝑡) and 𝐵(𝑡) are given by𝐴(𝑡)=𝐷_𝑥(𝑓)=[𝜕𝑓_𝑖𝜕𝑥_𝑗]∈𝑅^𝑛×𝑚,𝐵(𝑡)=𝐷_𝛼(𝑓)=[𝜕𝑓_𝑖𝜕𝛼_𝑗]∈𝑅^𝑛×𝑝. (A.2) Define the integrating factor 𝑀(𝑡,𝑠)=exp[^𝑡_𝑠𝐴(𝑢)𝑑𝑢]∈𝑅^𝑛×𝑛. (A.3) premultiplying both sides of (A.1) by 𝑀⁻¹(𝑡,0) and integrating, we get the solution of the linear nonhomogeneous and nonautonomous equation (A.1) as 𝛿𝑥(𝑡)=𝑀(𝑡,0)𝛿𝑥(0)+𝑀(𝑡,0)𝐿(𝑡,0)𝛿𝛼, (A.4) where 𝐿(𝑡,𝑠)=^𝑡_𝑠𝑀⁻¹(𝑢,0)𝐵(𝑢)𝑑𝑢∈𝑅^𝑛×𝑝. (A.5) From definitions (A.3)–(A.5), it can be verified that, for 𝑢<𝑠<𝑡,𝑀(𝑡,𝑢)=𝑀(𝑡,𝑠)𝑀(𝑠,𝑢),𝐿(𝑡,𝑢)=𝐿(𝑡,𝑠)+𝐿(𝑠,𝑢). (A.6) We now consider two special cases.

Case A. Let 𝛿𝛼=0, that is, the initial perturbations are confined only to the initial condition, 𝑥(0). Then setting 𝐵(𝑡)≡0, from (A.5) we see that 𝐿(𝑡,0)≡0. From (A.4) we get 𝛿𝑥(𝑡)=𝑀(𝑡,0)𝛿𝑥(0). (A.7)

Case B. Let 𝛿𝑥(0)=0, that is, the initial perturbations are confined only to the parameter, 𝛼. Then setting 𝐴(𝑡)≡0, from (A.3), we see that 𝑀(𝑡,𝑠)≡𝐼, the identity matrix. Now from (A.4) it follows that 𝛿𝑥(𝑡)=𝐿(𝑡,0)𝛿𝛼. (A.8)

B. Computation of Sensitivity Functions

Let 𝑐=(𝑐₁,𝑐₂,…,𝑐_𝑞)^𝑇∈𝑅^𝑞, where 𝑥=𝑥(𝑡,𝑐)=(𝑥₁(𝑡,𝑐),𝑥₂(𝑡,𝑐),…,𝑥_𝑛(𝑡,𝑐))^𝑇∈𝑅^𝑛 and ℎ(𝑥)=(ℎ₁(𝑥),ℎ₂(𝑥),…,ℎ_𝑚(𝑥))^𝑇∈𝑅^𝑚. Let 𝑎=(𝑎₁,𝑎₂,…,𝑎_𝑚)^𝑇∈𝑅^𝑚 and consider𝜙₁(𝑐)=𝑎^𝑇ℎ(𝑥)=^𝑚_𝑖=1𝑎_𝑖ℎ_𝑖(𝑥). (B.9) By applying the standard chain rule it can be verified that the gradient ∇_𝑐𝜙₁(𝑐) with respect to 𝑐 is given by∇_𝑐𝜙₁(𝑐)=𝐷^𝑇_𝑐(𝑥)𝐷^𝑇_𝑥(ℎ)𝑎, (B.10) where 𝐷_𝑥(ℎ)=[𝜕ℎ_𝑖𝜕𝑥_𝑗]∈𝑅^𝑚×𝑛 (B.11) is the Jacobian of ℎ and𝐷_𝑐(𝑥)=[𝜕𝑥_𝑖𝜕𝑐_𝑗]∈𝑅^𝑛×𝑞 (B.12) is the (Jacobian) sensitivity of the vector 𝑥 with respect to 𝑐 at time 𝑡.

Now consider a quadratic form 𝜙₂(𝑐)=ℎ^𝑇(𝑥)𝐴ℎ(𝑥), (B.13) where 𝐴∈𝑅^𝑚×𝑚 is a symmetric and positive definite matrix.

Then by the product rule ∇_𝑐𝜙₂(𝑐)=∇_𝑐[ℎ^𝑇𝑏]+∇_𝑐[𝑎^𝑇ℎ], (B.14) where 𝑏=𝐴ℎ=𝑎. By applying (B.10) to each of the two terms on the right side of (B.14), it follows that ∇_𝑐𝜙₂(𝑐)=2𝐷^𝑇_𝑐(𝑥)𝐷^𝑇_𝑥(ℎ)𝐴ℎ. (B.15) Finally, if 𝑄(𝑐)=[𝑧−ℎ(𝑥)]^𝑇𝐴[𝑧−ℎ(𝑥)], (B.16) then expand𝑄(𝑐)=𝑧^𝑇𝐴𝑧−2𝑧^𝑇𝐴ℎ(𝑥)+ℎ^𝑇(𝑥)𝐴ℎ(𝑥) (B.17) since 𝑧 does not depend on 𝑐, by using (B.10) and (B.15), we get𝑄(𝑐)=−2𝐷^𝑇_𝑐(𝑥)𝐷^𝑇_𝑥(ℎ)𝐴[𝑧−ℎ(𝑥)]. (B.18) Define𝜉=𝐷^𝑇_𝑥(ℎ)𝐴[𝑧−ℎ(𝑥)]∈𝑅^𝑛, (B.19) then∇_𝑐𝑄(𝑐)=−2𝐷^𝑇_𝑐(𝑥)𝜉. (B.20) That is, the gradient of 𝑄 with respect to 𝑐 is the linear combination of the sensitivity vectors, that is, the columns of 𝐷^𝑇_𝑐(𝑥), where the coefficients in this linear combination are the elements of the vector 𝜉.

C. Computation of ∇_𝛼𝐽 in (27)

Given 𝜂(𝑡_𝑖), 𝑀^𝑇(𝑡_𝑖,𝑡_𝑖−1),𝐿^𝑇(𝑡_𝑖,𝑡_𝑖−1) for 1≤𝑖≤𝑁, the expression on the right-hand side of (27) can be efficiently computed as shown in Algorithm 1.

DO 𝑗 = 𝑁 to 1  DO 𝑖 = 𝑁 to 𝑗  𝜂(𝑡𝑖)=𝑀𝑇(𝑡𝑗,𝑡𝑗−1)𝜂(𝑡𝑖)  END END 𝜇(1)=𝐿(𝑡1,0) DO 𝑖 = 2 to 𝑁 𝜇(𝑖)=𝜇(𝑖−1)+𝐿(𝑡𝑖,𝑡𝑖−1) END Grad = 0 DO 𝑖 = 1 to 𝑁 Grad = Grad + 𝜇(𝑖)𝜂(𝑡𝑖) END

Algorithm 1 

Then ∇_𝛼𝐽 = − Grad. It is to be noticed that there is only matrix-vector multiplication involved in these operations and not matrix-matrix multiplication.

Acknowledgments

At an early stage of formulating our ideas on this method of data assimilation, Fedor Mesinger and Qin Xu offered advice and encouragement. Qin Xu and Jim Purser carefully checked the mathematical development, and suggestions from the following reviewers went far to improve the presentation: Yoshi Sasaki, Bill Stockwell, and anonymous formal reviewers of the manuscript. S. Lakshmivarahan’s work is supported in part by NSF EPSCoR RII Track 2 Grant 105-155900 and by NSF Grant 105-15400, and J. Lewis acknowledges the Office of Naval Research (ONR), Grant No. N00014-08-1-0451, for research support on this project.