A data space search is the most efficient way to solve a linear, least-squares smoothing problem defined over a fixed time interval. The method exploits linearity, and so is unavailable for nonlinear dynamics, or for penalties other than least-squares. As discussed in Chapter 3, a data space search may be conducted on linear iterates of the nonlinear Euler–Lagrange equations. The existence of the nonlinear equations implies that the penalty is a smooth functional of the state, in which case a state space search may always be initiated. The nature of state space searches is intuitively clear, and their use is widespread. Conditioning degrades as the size of the state space gets very large. Collapsing the size of the state space by assuming “perfect” dynamics is the basis of “the” variational adjoint method: only initial values, boundary values and parameter values are varied. Preconditioning may in principle be effected by use of second-order variational equations, but even iterative construction of the state space preconditioner is unfeasible for highly realistic problems. Technique for numerical integration of variational equations is not a paramount consideration, but deserving of attention since it can be consuming of human time.

Operational forecasting is inherently sequential; data are constantly arriving and forecasts must be issued regularly. In such an environment, it is more natural to filter a model and data sequentially than to smooth them over a fixed interval. […]

Chapter 1 is a minimal course on assimilating data into models using the calculus of variations. The theory is introduced with a “toy” model in the form of a single linear partial differential equation of first order. The independent variables are a spatial coordinate, and time. The well-posedness of the mixed initial-boundary value problem or “forward model” is established, and the solution is expressed explicitly with the Green's function. The introduction of additional data renders the problem ill-posed. This difficulty is resolved by seeking a weighted least-squares best fit to all the information. The fitting criterion is a penalty functional that is quadratic in all the misfits to the various pieces of information, integrated over space and time as appropriate. The best-fit or “generalized inverse” is expressed explicitly with the representers for the penalty functional, and with the Green's function for the forward model. The behavior of the generalized inverse is examined for various limiting choices of weights. The smoothness of the inverse is seen to depend upon the nature of the weights, which will be subsequently identified as kernel inverses of error covariances. After reading Chapter 1, it is possible to carry out the first four computing exercises in Appendix A.

It is a long road from deriving the formulae for the generalized inverse of a model and data to seeing results. First experiments (McIntosh and Bennett, 1984) involved a linear barotropic model separated in time, simple coarsely-resolved numerical approximations, a handful of pointwise measurements of sea level and a serial computer. Contemporary models of oceanic and atmospheric circulation involve nonlinear dynamics and parameterizations, advanced high-resolution numerical approximations, vast quantities of data often of a complex nature, and parallel computers. Chapter 3 introduces some general principles for travelling this long road of implementation.

The first principle is accelerating the representer algorithm by task decomposition, that is, by simultaneous computation of representers on parallel processors. The objective may be either the full representer matrix as required by the direct algorithm, or a partial matrix for preconditioning the indirect algorithm. The calculation of an individual representer, or indeed any backward or forward integration, may itself be accelerated by domain decomposition, but this is a common challenge in modern numerical computation (Chandra et al., 2001; Pacheco, 1996) and will not be addressed here. Even without considering the coarse grain of task decomposition or the fine grain of domain decomposition, the direct and indirect representer algorithms for linear inverses are highly intricate. Schematics are provided here in the form of “time charts”.

Dynamical errors and input errors may be correlated in space or in time or in both. Error covariances must be convolved with adjoint variables. […]

An ocean data assimilation system in miniature

The pages of this book are filled with the mathematics of oceanic and atmospheric circulation models, observing systems and variational calculus. It would only be natural to ask: What is going on here, and is it really new? The answers are “regression” and hence “no”: almost every issue of any marine biology journal contains a variational ocean data assimilation system in miniature.

Linear regression in marine biology

The article “Repression of fecundity in the neritic copepod Acartia clausi exposed to the toxic dinoflagellate Alexandrium lusitanicum: relationship between feeding and egg production”, by Jörg Dutz, appeared in Marine Ecology Progress Series in 1998.

The calculus of variations uses Green's functions and representers to express the best fit to a linear model and data. Mathematical construction of the representers is devious, and the meaning of the representer solution to the “control problem” of Chapter 1 is not obvious. There is a geometrical interpretation, in terms of observable and unobservable degrees of freedom. Unobservability defines an orthogonality, and the representers span a finite-dimensional subspace of the space of all model solutions or “circulations”. The representers are in fact the observable degrees of freedom.

A statistical interpretation is also available: if the unknown errors in the model are regarded as random fields having prescribed means and covariances, then the representers are related, via the measurement processes, to the covariances of the circulations. Thus the representer solution to the variational problem is also the optimal linear interpolation, in time and space, of data from multivariate, inhomogeneous and nonstationary random fields. The minimal value of the penalty functional that defines the generalized inverse or control problem is a random number. It is the χ2 variable, if the prescribed error means or covariances are correct, and has one degree of freedom per datum. Measurements need not be pointwise values of the circulation; representers along with their geometrical and statistical interpretations may be constructed for all bounded linear measurement functionals.

Analysis of the conditioning of the determination of the representer amplitudes reveals those degrees of freedom which are the most stable with respect to the observations. […]

The “toy” forward model introduced in Chapter 1 defines a well-posed mixed initial value–boundary value problem. The associated operator (wave operator plus initial operator plus boundary operator) is invertible, or nonsingular. Specifying additional data in the interior of the model domain renders the problem overdetermined. The operator becomes uninvertible, or singular. The difficulty may be resolved by constructing the generalized inverse of the operator, in a weighted least-squares sense.

An important class of regional models of the ocean or atmosphere defines an ill-posed initial-boundary value problem, regardless of the choice of open boundary conditions. All flow variables may as well be specified on the open boundaries. The excess of information may be regarded as data on a bounding curve, rather than at an interior point. The difficulty may again be resolved by constructing the generalized inverse in the weighted least-squares sense. The Euler–Lagrange equations form a well-posed boundary value problem in space–time. Solving them by forward and backward integrations is precluded, since no partitioning of the variational boundary conditions yields well-posed integrations. The penalty functional must be minimized directly.

Things are different if the open region is moving with the flow.

These computational exercises complement the analytical development of variational data assimilation in the text, and also serve to develop the confidence needed for more ambitious calculations. All code for these exercises is available at an anonymous ftp site. The linear, one-dimensional “toy” model of §1.1 is upgraded here to a linear, two-dimensional shallow-water model in both continuous and finite-difference form. Continuous and discrete penalty functionals are developed, and the respective Euler–Langrange equations are derived. Representers are calculated directly, so the generalized inverse may then be calculated directly or indirectly.

Inverse modeling has many applications in oceanography and meteorology. Charts or “analyses” of temperature, pressure, currents, winds and the like are needed for operations and research. The analyses should be based on all our knowledge of the ocean or atmosphere, including both timely observations and the general principles of geophysical fluid dynamics. Analyses may be needed for flow fields that have not been observed, but which are dynamically coupled to observed fields. The data must therefore contribute not only to the analyses of observed fields, but also to the inference of corrections to the dynamical inhomogeneities which determine the coupled fields. These inhomogeneities or inputs are: the forcing, initial values and boundary values, all of which are themselves the products of imperfect interpolations. In addition to input errors, the dynamics will inevitably contain errors owing to misrepresentations of phenomena that cannot be resolved computationally; the data are therefore also required to improve the dynamics by adjusting the empirical coefficients in the parameterizations of the unresolved phenomena. Conversely, the model dynamics must have some credibility, and should be allowed to influence assessments of the effectiveness of observing systems. Finally, and perhaps most compelling of all, geophysical fluid dynamical models need to be formulated and tested as formal scientific hypotheses, so that the development of increasingly realistic models may proceed in an orderly and objective fashion. All of these needs can be met by inverse modeling.

A brief history of bounds derived using the analytic method

Bergman (1978) recognized that the analytic properties discussed in chapter 18 on page 369 provide a powerful tool for deriving bounds. He rederived the Hashin-Shtrikman bounds and obtained new bounds correlating different properties of composites. A major success of the approach was that it lead to tight bounds on the complex dielectric constant of a two-phase composite (Milton 1979, 1980, 1981a; Bergman 1980, 1982). These bounds are illustrated in figure 27.1 on the next page. [The first available bounds on complex dielectric constants were those of Schulgasser and Hashin (1976), but they were limited to materials with low-loss constituents, that is, with permittivities having small imaginary parts.] These complex dielectric constant bounds have been directly compared with experimental measurements: Niklasson and Granqvist (1984) applied them to bounding the optical properties of composite films; Korringa and LaTorraca (1986) applied them to bounding the complex electrical permittivity of rocks; Golden (1995) applied them to bounding the complex permittivity of sea ice; and Mantese, Micheli, Dungan, Geyer, Baker-Jarvis, and Grosvenor (1996) applied them to bounding the complex dielectric constant and magnetic permeability of composites of Barium Titanate and ferrite. In most cases the experimental measurements were consistent with the bounds. However, it is important to recognize that these bounds apply only in the quasistatic limit where the wavelength of the radiation is much larger than the inhomogeneities of the microstructure; see Aspnes (1982). McPhedran, McKenzie, and Milton (1982); McPhedran and Milton (1990); and Cherkaeva and Golden (1998) applied the bounds in an inverse fashion to obtain quite tight bounds on the volume fraction from measurements of the complex dielectric constant.

Variational principles have long been known in the context of both conductivity type problems and elasticity problems. Their application to composites was initiated by Hill (1952), who used them to show that the Voigt (1889, 1910) and Reuss (1929) estimates of the elastic moduli of polycrystals were in fact bounds. A new type of variational principle was discovered by Hashin and Shtrikman (1962a, 1962b, 1963), which become famous because it lead to optimal bounds on the conductivity, bulk, and shear moduli of isotropic composites of two isotropic phases. Hill (1963b) gave rigorous proof of their variational principles, showing how they could be derived from the classical variational principles. Subsequently, Hashin (1967) generalized the variational principles to inhomogeneous elastic bodies (not just composites) subject to body forces and mixed boundary conditions. Cherkaev and Gibiansky (1994) extended all of these variational principles to media with complex moduli. Talbot and Willis (1985) extended the Hashin and Shtrikman variational principles to nonlinear media. Other variational inequalities for nonlinear media, based on comparisons with linear inhomogeneous media, were obtained by Ponte Castañeda (1991). In a development that falls outside the range of this book, Smyshlyaev and Fleck (1994, 1996) extended the Hashin and Shtrikman variational principles to elastic composites, where the elastic energy depends not only on the strain, but also on the strain gradient.

Classical variational principles and inequalities

We have seen how to manipulate equations into a form where the tensor entering the constitutive law is self-adjoint and positive-definite.