Abstract

A method is presented in which the adjoint of a four-dimensional variational data assimilation system (4D-Var) was used to compute the expected analysis and forecast error variances of linear functions of the ocean state vector. The power and utility of the approach are demonstrated using the Regional Ocean Modeling System configured for the California Current system. Linear functions of the ocean state vector were considered in the form of indices that characterize various aspects of the coastal upwelling circulation. It was found that for configurations of 4D-Var typically used in ocean models, reliable estimates of the expected analysis error variances can be obtained both for variables that are observed and unobserved. In addition, the contribution of uncertainties in the model control variables to the forecast error variance was also quantified. One particularly powerful and illuminating aspect of the adjoint 4D-Var approach to the forecast problem is that the contribution of individual observations to the predictability of the circulation can be readily computed. An important finding of the work presented here is that despite the plethora of available satellite observations, the relatively modest fraction of in situ subsurface observations sometimes exerts a significant influence on the predictability of the coastal ocean. Independent checks of the analysis and forecast error variances are also presented, which provide a direct test of the hypotheses that underpin the prior error and observation error estimates used during 4D-Var.

1. Introduction

Data assimilation forms a critical component of all modern analysis and forecast systems for both the atmosphere and ocean. As such, an important component of all data assimilation systems should be an estimate of the posterior or analysis error covariance associated with the best circulation estimate, as well as an estimate of the ensuing forecast error covariance. In addition, when data assimilation is applied sequentially, the resulting analysis typically becomes the prior for the next analysis cycle. Therefore, in addition to providing information about the expected uncertainty in the analysis, the analysis error covariance can be used to update the prior or background error covariance for the next data assimilation cycle.

While analysis error information can in principle be readily computed for statistical interpolation schemes and Kalman filters (Daley 1991), it presents a challenge for variational methods because the information needed is generally not available in a convenient form. In this paper, we will present a method by which formal posterior covariance information can be obtained using the adjoint of the entire variational data assimilation algorithm. To set the stage for our work, we will begin with a brief review of four-dimensional variational assimilation (4D-Var) of ocean data; however, the same ideas obviously apply to the atmosphere as well, and to 3D-Var.

Following the notation introduced by Ide et al. (1997) and expanded by Daget et al. (2009) and Moore et al. (2011a, hereafter M11a), we denote by x the state vector of the ocean (or atmosphere). Data assimilation seeks to combine, in an optimal way, a background (or prior) estimate of the circulation, xb(t), with observations that are arranged in the vector yo. The background circulation is generally a solution of a numerical model, and depends on background estimates of the initial conditions, xb(t0); surface forcing, fb(t); and open boundary conditions, bb(t). The model solution will be denoted by xb(ti) = M(ti, ti−1)[xb(ti−1), fb(ti), bb(ti)], where M(ti, ti−1) represents the model operators and fb(ti) and bb(ti) represent the actions of the forcing and boundary conditions, respectively, during the interval [ti−1, ti]. During 4D-Var the best circulation estimate is identified by adjusting a vector of control variables composed of the initial conditions, surface forcing, open boundary conditions, and corrections for model error. Sometimes, 4D-Var is implemented using the incremental formulation of Courtier et al. (1994), in which the control vector, denoted δz, is composed of increments to the priors, specifically increments δx(t0) to the state vector at initial time t0, as well as increments to surface forcing, δf(t); open boundary conditions, δb(t); and corrections for model error, η(t) at all times t spanning the data assimilation cycle. To first order, the state vector increments δx(t) are governed by δx(ti) = (ti, ti−1)δu(ti−1), where δu(ti−1) = [δxT(ti−1), δfT(ti), δbT(ti), ηT(ti)]T, and (ti, ti−1) represents the perturbation tangent linear model linearized about the time-evolving background xb(t).

Denoting as zb = [xbT(t0), …, fbT(ti), …, …, bbT(ti), …, …, ηbT(ti), …]T the vector of background fields (where the ellipses denote the sequence of all times ti in the assimilation interval), then the control vector that yields the best circulation estimate can be conveniently expressed as , where δza = [δxT(t0), …, δfT(ti), …, …, δbT(ti), …, …, δηT(ti), …]T is the vector of analysis (or posterior) increments that minimize the cost function:

 
formula

where and are the background and observation error covariance matrices, respectively; d = {yoH[xb(t)]} is the innovation vector; H[xb(t)] is the background evaluated at the observation points via the observation operator H; and represents the tangent linear model sampled at the observation points, and is essentially a convolution in time of (ti, ti−1) and the tangent linearization of H.

The analysis increments can be expressed as δza = d, where

 
formula
 
formula

is the gain matrix. Equation (2) is sometimes referred to as the primal form of and arises from searching for δza in the full space spanned by the control vector δz. Alternatively, (3) is referred to as the dual form of and arises from a search for δza in the space of linear functions of δz spanned only by the observations (i.e., δz), usually referred to as observation space.

Using (2) or (3), the posterior or analysis error covariance matrix, a, of the best estimate can be written in several different (although equivalent) forms; namely,

 
formula
 
formula
 
formula

It is important to realize that the analysis error covariance matrix given by (4)(6) will only be the true a if the prior error covariance matrix and the observation error covariance matrix are correctly specified, and if the assumptions underlying the analysis system are correct, namely Gaussian, unbiased errors in the priors and the observations.

Equation (4) shows that a is equal to the inverse of the Hessian matrix, , of the cost function J in (1). During applications of primal 4D-Var, the leading eigenvalues and eigenvectors of can be estimated and used to compute a reduced-rank approximation of a. However, since a = −1, the leading eigenvectors of are the least important eigenvectors of a and, in practice, explain very little of the analysis error variance. Nonetheless, while such estimates may be expected to overestimate the analysis error, they are used routinely (e.g., Fisher and Courtier 1995) and provide information about geographic variations of the error. Equations (5) and (6) arise directly from the dual form of in (3), and Daley and Barker (2001) [see also Gelaro et al. (2002) for more details] describe a routine approach for estimating a based on a Cholesky factorization of the stabilized representer matrix .

In practice, the matrix inverse in (2) or (3) is evaluated iteratively by solving an equivalent system of linear equations using a conjugate gradient algorithm. In a number of 4D-Var systems in current use (e.g., Fisher and Courtier 1995; Tshimanga et al. 2008; Chua et al. 2009; M11a), the conjugate gradient algorithm is formulated in terms of a sequence of Lanczos vectors (Lanczos 1950), which are used to construct a reduced-rank approximation of . By virtue of the orthonormal property of the Lanczos vectors, (6) leads to a convenient reduced-rank form for a. Moore et al. (2011b, hereafter M11b) have used (6) in this capacity in an oceanographic application of 4D-Var to estimate analysis errors for the California Current circulation. However, since the Lanczos basis employed does not typically span the full space of , the resulting expected analysis error is still overestimated.

Monte Carlo methods (Bennett 2002) and ensemble techniques (Houtekamer et al. 1996; Žagar et al. 2005; Belo Pereira and Berre 2006; Berre et al. 2006) can also be used to provide estimates of a in conjunction with 4D-Var, and the method we describe in this paper is inspired by this idea. In section 2 we describe how the adjoint of a 4D-Var system can be used to estimate the expected analysis error covariance without the need to explicitly compute an ensemble. The use of the adjoint of a data assimilation system is not new, and has traditionally been used in numerical weather prediction to provide observation sensitivity and observation impact information pertaining to analyses (e.g., Zhu and Gelaro 2008) and forecasts (e.g., Langland and Baker 2004; Gelaro and Zhu 2009). In this study, however, the utility of the adjoint of a data assimilation system has been extended to provide formal error variance information about theoretical infinite-sized ensembles. Both our approach and those of previous investigations, however, capitalize on the property of the adjoint of the data assimilation system to provide information about the sensitivity of scalar functions of the control vector to variations in the inputs to the system, in this case the observations and priors. The method has been tested in relation to the circulation of the California Current system (CCS), and the configuration of the model and 4D-Var system are presented in section 3. The adjoint of 4D-Var is used to estimate the expected errors associated with various dynamical aspects of the coastal upwelling circulation that develops along the central California coast, and indices that quantify this circulation are introduced in section 4. The expected analysis and forecast error variance of the circulation indices are explored in sections 5 and 6, respectively. Section 7 explores the predictability of the circulation attributable to the different observation platforms, and in section 8 we present a series of checks for consistency and statistical reliability for the derived error estimates. We end with a summary and conclusions in section 9.

2. Analysis error covariance from perturbed analyses

Consider an ensemble of analyses that are generated by perturbing the observations and background fields with perturbations drawn from normal distributions with covariances and , respectively. It has been proven (e.g., Žagar et al. 2005; Belo Pereira and Berre 2006; Daget et al. 2009) that the covariance of the perturbed analyses about the unperturbed analysis simulates the analysis error covariance. If the prior covariances and are equal to the true covariances, the covariance of the analysis ensemble will equal the true analysis error covariance. This is the principle behind ensemble 4D-Var currently used operationally at some numerical weather prediction (NWP) centers (e.g., Isaksen et al. 2010), and is the underlying premise behind the method presented here. In our case, however, instead of explicitly generating an ensemble of analyses, we use the adjoint of the analysis system in conjunction with the expectation operator and consider all possible analyses drawn from the appropriate distribution.

As shown by M11a (see also the  appendix), the circulation analysis can be alternatively expressed as , where the analysis increment δza = (d) and (d) represents the entire 4D-Var procedure, which, in general, is a nonlinear function of the innovation vector d by virtue of the algorithm employed to minimize J in (1). The ensemble 4D-Var approach outlined above involves perturbing the observations and background control vector, which yields perturbations δd in the innovations, and associated first-order changes in the analysis increments given by (∂/∂d)δd, where ∂/∂d represents the tangent linearization of the entire 4D-Var procedure. With this in mind, consider a 4D-Var analysis spanning the interval [t0, t0 + τ]. The time evolution of the resulting posterior state vector is given by

 
formula

Consider now an ensemble of 4D-Var analyses , where the lth ensemble member is generated by adding perturbations to the background control vector zb and perturbations to the observations, where and are drawn from normal distributions with covariances and , respectively. The resulting perturbation in the 4D-Var analysis control vector is given by , where is the perturbation to the innovation vector. The tangent linear 4D-Var operator ∂/∂d is linearized about the unperturbed 4D-Var sequence of iterations that yield the unperturbed analysis xa, fa, ba, and ηa. We can also express the analysis perturbation control vector as , so that from (7) the state vector of each perturbed analysis evolves according to

 
formula

where a is the tangent linear model linearized about the unperturbed time-evolving analysis (7). Since our primary focus in later sections is on the state vector x, we introduce an alternative form of the tangent linear model operator, a, that isolates the state vector perturbation δx arising from a control vector perturbation δz, so that δx(ti) = a(ti, t0)δz. Recall that is the vector of perturbations in the analysis at a single instant in time ti, while δzl are the perturbations in all the control vector elements (i.e., δx at t0 and δf, δb, and δη at all times during the assimilation interval). At the end of the analysis cycle, the perturbation in the analysis state vector will to first order be given by

 
formula

As noted earlier, if 4D-Var is run to complete convergence, the priors and are correctly specified, and all of the underlying assumptions about 4D-Var are correct (e.g., an unbiased estimator), then the variance of the perturbed 4D-Var analyses (8) about the unperturbed analysis (7) equals the expected analysis error covariance. Therefore, using (9), the analysis error covariance of the state vector at time t0 + τ is given by

 
formula

where 〈⋯〉 denotes the expectation operator and (∂/∂d)T represents the adjoint of the tangent linearized 4D-Var algorithm. Formally, the tangent linear model a is linearized about the unperturbed 4D-Var analysis (indicated by the subscript a). However, in the case where the linearized dynamics are always used to propagate information, as in the case of the representer method of Bennett (1992) used in this study, then a b and the tangent linear model is linearized about the unperturbed background (indicated by the subscript b).

Computing the entire matrix given by (10) presents a considerable challenge due to its large dimension. However, for linear scalar functions = [x(t)] (e.g., space–time integrals of x), computation of the expected analysis error variance of presents much less of a challenge. Consider the generic linear scalar function of the discrete model state vector xkx(t0 + kΔt) given by , where Δt is the model time step and hk is an appropriate vector. The prior error variance of (i.e., before data assimilation) is given by

 
formula

where (x)k x(t0 + kΔt, t0) is given by a or b depending on the 4D-Var algorithm used (see below). Similarly, the expected posterior error variance in after data assimilation is given by

 
formula

where . If 4D-Var is run to complete convergence, then for a linear data assimilation system (or the case of one outer loop; see section 4) ∂/∂d = , and is identical to (5) for the true gain matrix.

The model used here is the Regional Ocean Modeling System (ROMS), for which an extensive suite of 4D-Var tools has recently been developed (M11a) based on both the primal and dual formulations. Two dual 4D-Var algorithms are currently available in ROMS: one that is based on the Physical-space Statistical Analysis System (PSAS; Cohn et al. 1998), and one that utilizes the indirect representer method of Egbert et al. (1994), referred to here as R4D-Var. In the case of primal 4D-Var or dual PSAS, the posterior circulation estimate is a solution of the nonlinear model, in which case xa in (11) and (12). In the case of the dual indirect representer method, the posterior circulation estimate is a solution of a finite-amplitude linearization of the model (Bennett 2002) linearized about the prior circulation, in which case xb in (11) and (12). In the present study we will focus our attention on R4D-Var.

3. Configuration of ROMS and 4D-Var

ROMS is a hydrostatic, primitive-equation, Boussinesq ocean general circulation model designed primarily for coastal applications by employing terrain-following vertical coordinates and horizontal orthogonal curvilinear coordinates (Shchepetkin and McWilliams 2005). The configuration of ROMS and R4D-Var used in the calculations reported below is described in detail by Broquet et al. (2009a,b, 2011) and M11b, so only a brief description will be given here.

The ROMS CCS domain spans the region 31°–48°N to 134°–116°W with 30-km horizontal resolution and 30 terrain-following σ levels in the vertical. While this resolution is marginal for modeling the mesoscale circulation of the CCS, it provides a reasonable representation of the large-scale circulation. The motivation for using a coarse resolution is that the sequential 4D-Var computations presented here are computationally very demanding. Since the purpose of this paper is to demonstrate the efficacy of using the adjoint of 4D-Var for estimating posterior errors of a sequential data assimilation system, the use of a coarse-resolution model is justified. The model domain and bathymetry are shown in Fig. 1.

Fig. 1.

The model domain and bathymetry. Also shown are the sections along which alongshore and cross-shore transports were computed, and the region in which upper-ocean heat content and SST were evaluated.

Fig. 1.

The model domain and bathymetry. Also shown are the sections along which alongshore and cross-shore transports were computed, and the region in which upper-ocean heat content and SST were evaluated.

The model forcing was derived from daily averaged output of atmospheric boundary layer fields from the Naval Research Laboratory’s (NRL) Coupled Ocean–Atmosphere Mesoscale Prediction System (COAMPS; Doyle et al. 2009) using the bulk formulations of Liu et al. (1979) and Fairall et al. (1996a,b). The resulting surface fluxes of momentum, heat, and freshwater represent the background surface forcing, denoted fb(t) in section 1. The model domain has open boundaries at the northern, southern, and western edges where the tracer and velocity fields were prescribed, and the free surface and vertically integrated flow were subject to Chapman (1985) and Flather (1976) boundary conditions, respectively. The prescribed open boundary solution was taken from the Estimating the Circulation and Climate of the Ocean (ECCO) global data assimilation product (Wunsch and Heimbach 2007) and represents the background open boundary conditions denoted bb(t) in section 1. A sponge layer was also used adjacent to each open boundary in which viscosity increased linearly from 4 m2 s−1 in the interior to 100 m2 s−1 at the boundary over a distance of ~100 km.

Following M11b, ROMS R4D-Var was run sequentially for the period July 2002–December 2004 starting from a set of background initial conditions on 27 July 2002 derived from the 4D-Var sequence of analyses computed by Broquet et al. (2009a). The observations assimilated in the model were collected by various platforms, and include gridded sea surface height (SSH) analyses in the form of dynamic topography from Aviso at ~⅓° resolution every 7 days; a blended SST product with 10-km resolution, available daily, and consisting of 5-day means derived from the Goddard Earth Observing System (GEOS), Advanced Very High Resolution Radiometer (AVHRR), and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite instruments [courtesy of D. Foley at the National Oceanic and Atmospheric Administration (NOAA) and available east of 130°W only at the time the experiments presented here were performed]; in situ hydrographic observations extracted from the quality controlled EN3 (version v1d) data archive maintained by the Met Office as part of the European Union ENSEMBLES project (Ingleby and Huddleston 2007); and tagged elephant seal data from the Tagging of Pacific Pelagics (TOPP) program. To reduce data redundancy, all observations within each model grid cell, over a 6-h time window, were combined to form “superobservations,” and the standard deviation of the observations that contribute to the superobservation in each grid cell was used as an estimate of the error of representativeness.

Observation errors were assumed to be uncorrelated in space and time, resulting in a diagonal observation error covariance matrix . The variances along the main diagonal of were assigned as a combination of measurement error and the error of representativeness. Measurement errors were chosen independent of the data source, with the following standard deviations: 0.02 m for SSH, 0.4°C for SST, 0.1°C for in situ T, and 0.01 for in situ S.

A background error standard deviation associated with xb(t0), fb(t), and bb(t) was estimated for each calendar month. Following Weaver et al. (2005), the background error for the initial conditions was decomposed into a balanced and unbalanced component, and the associated background error covariance matrix was factorized as the product of a univariate covariance matrix (describing the errors of the unbalanced flow) and a dynamical balance operator. The dynamical balance operator yields cross-covariance information about the errors and is based on the TS properties of the water column, hydrostatic balance, and geostrophic balance (Weaver et al. 2005; M11a). The background error standard deviations for the unbalanced components of the control vector were estimated each month based on the variance of the model run during the period 1999–2004 subject only to surface forcing (i.e., no data assimilation). The temporal variability of the COAMPS surface forcing for the period 1999–2004 was used as the variance for the background surface forcing error, and the open boundary condition background error variances were chosen to be the variances of the ECCO fields at the boundaries. A balance operator was not applied to the surface forcing or open boundary condition increments. In all of the experiments presented here, the circulation estimates were computed subject to the strong constraint, in which case the model is assumed to be free of error, and the corrections η for model error were set to zero.

At the present time, ROMS 4D-Var supports only homogeneous error correlations that are separable in the horizontal and vertical. The horizontal and vertical correlation functions are modeled using the diffusion operator approach described by Weaver and Courtier (2001) (see also M11a for details of the implementation in ROMS). The decorrelation length scales used to model the background errors of all initial condition control variable components of were 50 km in the horizontal and 30 m in the vertical. Horizontal correlation scales chosen for the background surface forcing error components of were 300 km for wind stress and 100 km for heat and freshwater fluxes. The correlation lengths for the background open boundary condition error components of were 100 km in the horizontal and 30 m in the vertical. Explicit temporal correlations of the background errors are not included in the current version of ROMS 4D-Var. However, the surface forcing and boundary condition increments, δf(t) and δb(t), were only computed daily and interpolated to each intervening model time step, which effectively introduces some temporal correlation of the errors in fb(t) and bb(t). A discussion of the choice of the aforementioned background error covariance parameters can be found in Broquet et al. (2009a,b, 2011) and M11b.

The R4D-Var model was run sequentially for the period July 2002–December 2004 with 1 outer loop and 60 inner loops1 spanning 7-day assimilation cycles, which is sufficient to guarantee a significant level of convergence of J toward its minimum value as demonstrated by M11b. Following El Akkraoui and Gauthier (2010), the minimum residual method of Paige and Saunders (1975) was used to minimize the cost function J during each data assimilation cycle, and the analysis at the end of each assimilation cycle, xa(t0 + 7), becomes the background initial condition for the next assimilation cycle. Figure 2 shows the ratio of the final and initial values of J for the R4D-Var sequence, and indicates that an order of magnitude reduction in J during each cycle is typical.

Fig. 2.

A time series of the ratio of the final value of the cost function, Jf, to the initial value of the cost function, Ji, for the sequence of R4D-Var data assimilation cycles spanning the period July 2002–December 2004.

Fig. 2.

A time series of the ratio of the final value of the cost function, Jf, to the initial value of the cost function, Ji, for the sequence of R4D-Var data assimilation cycles spanning the period July 2002–December 2004.

4. Coastal transport, heat content, and SST

The central California coastal circulation is characterized by a pronounced seasonal cycle of coastal upwelling (Checkley and Barth 2009). Equatorward of Cape Mendocino (about 40°N), the winds are alongshore and upwelling favorable year round. However, during late spring and early summer as the atmospheric subtropical high pressure system builds over the ocean, the alongshore winds intensify and the coastal upwelling reaches a peak. With these characteristics of the circulation in mind, we will consider four linear functions (x). These functions also serve to demonstrate how (10) can provide error information not only about ocean-state variables that are frequently and directly observed (such as ocean temperature) but also error information about variables that are not directly observed. The following four functions are considered.

a. Alongshore transport crossing 37°N, 37N

The pronounced seasonal changes in the alongshore winds near the central California coast modulate the strength of the California Current system, as well as the cross-shore pressure gradients that develop in response to Ekman divergence at the coast and Ekman pumping associated with wind stress curl. As a representative measure of these changes, we consider the 7-day-average transport crossing 37°N between the coast and 127°W over the upper 500 m of the water column during each R4D-Var assimilation cycle. This section is indicated in Fig. 1. Specifically, following the notation introduced in section 2, we will denote this transport over the interval [t0, t0 + 7] as , where t0k1Δt and t0 + 7 ≡ k2Δt. The elements of the vector (h37N)k are zero except for those that correspond to the velocity grid points that contribute to the transport normal to the 37°N section for the prescribed longitude and depth range.

Figure 3a shows a time series of 37N(xb) computed from the background circulation xb(t) during each R4D-Var cycle. A pronounced seasonal cycle is apparent with poleward flow during winter and spring, and equatorward flow during summer and fall, and a peak-to-peak transport of about ±5 Sv (1 Sv ≡ 106 m3 s−1). Also shown in Fig. 3a is the increment in 37°N transport Δ37N = 37N(xa) − 37N(xb), the transport difference between the analysis, xa(t) and the background xb(t) due to assimilating the observations during each R4D-Var cycle. The mean transport increment is −0.35 SV with a range of ~±2 SV.

Fig. 3.

Time series of (a) 37N (solid line) and Δ37N (line with circles), (b) 500m (solid line) and Δ500m (line with circles), (c) HC (solid line) and ΔHC (line with circles), and (d) SST (solid line) and ΔSST (line with circles).

Fig. 3.

Time series of (a) 37N (solid line) and Δ37N (line with circles), (b) 500m (solid line) and Δ500m (line with circles), (c) HC (solid line) and ΔHC (line with circles), and (d) SST (solid line) and ΔSST (line with circles).

b. Cross-shore transport at the continental shelf break, 500m

The wind stress along the coast and wind stress curl farther offshore drive Ekman transport and Ekman pumping in the upper ocean, resulting in a cross-shore transport. As a representative measure of these processes during each R4D-Var assimilation cycle, we considered the 7-day-average transport in the upper 15 m of the water column, between 35° and 40.5°N, crossing over the continental shelf break, taken here to be the 500-m isobath. This section is indicated in Fig. 1. Specifically we will denote this transport over the interval [t0, t0 + 7] as , where the elements of the vector (h500m)k are zero except for those corresponding to the velocity grid points that contribute to the transport normal to the 500-m isobath for the prescribed longitude, latitude, and depth range. The vector (h500m)k is defined so that offshore transport corresponds to 500m > 0. Figure 3b shows time series of 500m(xb) computed from the background circulation xb(t) of each R4D-Var cycle, and shows that the transport is generally offshore during most of the year, with peak-to-peak variations ~±1 Sv. The cross-shore transport increments Δ500m = 500m(xa) − 500m(xb) are also shown in Fig. 3b, and have a mean close to zero, as well as peak-to-peak variations of ~±0.5 Sv.

c. SST and upper-ocean heat content, SST and HC

The seasonal cycle of coastal upwelling leads to pronounced changes in the thermal structure of the upper ocean. To quantify these changes, we consider the 7-day mean SST averaged over the target area shown in Fig. 1, and given by , where the elements of (hSST)k are zero except for those surface temperature grid points that fall within the target area. The target area encompasses the region where the largest seasonal variations in SST occur associated with the seasonal variations in the wind. In addition, we also consider the 7-day mean temperature in the same target area averaged over the upper 50 m of the water column; namely, , where (hHC)k is appropriately defined. The linear function HC is also proportional to the heat content in the same volume of water. Time series of HC(xb) and SST(xb) for the R4D-Var background circulation are shown in Figs. 3c and 3d, respectively. The two time series are similar, and show that the lowest temperatures associated with the peak in upwelling occur during the spring and early summer. Following the weakening of the subtropical high, the alongshore upwelling-favorable winds decrease in strength around midsummer and the temperature of the surface waters increases. The increments ΔHC and ΔSST are also shown in Figs. 3c and 3d, and the mean increment in both cases is ~0.3°C, indicating a cold bias in the background. While HC(xb) and SST(xb) provide similar information, important differences in their error variance properties will become apparent later.

5. Posterior error estimates

The prior and posterior error variance of each linear function can be computed using (11) and (12), and Fig. 4 shows time series of the expected error variances of the various circulation indices introduced in section 4. Figure 4 shows that for all , is consistently lower than , indicating the expected positive impact of data assimilation. Recall that will only be the true error variance if and are correctly specified. Therefore, compared to the prior estimate, it is possible that the posterior estimate might actually be degraded by assimilating data if the and are not correct, yet the expected may indicate otherwise. Nonetheless, in the case of the transports 37N and 500m (Figs. 4a and 4b), the posterior error is ~40%–80% lower than the prior error, and both and tend to be larger during winter, a reflection of the seasonal dependence imposed on (verified by comparison with cases where there were no seasonal variations in ; not shown). In the case of HC and SST, the are two orders of magnitude smaller than the .

Fig. 4.

Time series of the background (prior) error variance (solid line) and the expected analysis (posterior) error computed from (12) (line with closed circles), and computed using the practical gain matrix (line with open circles) for (a) 37N, (b) 500m, (c) HC, and (d) SST.

Fig. 4.

Time series of the background (prior) error variance (solid line) and the expected analysis (posterior) error computed from (12) (line with closed circles), and computed using the practical gain matrix (line with open circles) for (a) 37N, (b) 500m, (c) HC, and (d) SST.

Figure 4 also shows an estimate of the posterior error computed using the reduced-rank approximation of the gain matrix according to , where . As discussed in the  appendix, can be computed using the Lanczos vectors derived from the m inner loops employed in R4D-Var. However, Fig. 4 indicates that is only marginally less than for the transports, and is significantly larger than in the case of HC and SST. This is because is based on mNobs, the number of observations, and represents a single realization of that typically will span only a small subspace of observation space. As such, tends to overestimate the true expected posterior error variance as discussed in M11b. However, even though the tangent linearization of 4D-Var, (∂/∂d), is computed using the same number of inner loops mNobs, each realization of the perturbation in the analysis increment (∂/∂d)δdl will span a different set of m directions in observation space. In the case of (12), all possible realizations of δdl, drawn from a normal distribution with covariance , are considered by virtue of the expectation operator used to derive , and represents a more reliable estimate of the true expected posterior error variance than , as is demonstrated in the  appendix.

6. Forecast errors

Consider now the adjacent analysis and forecast cycles illustrated schematically in Fig. 5a, where forecast cycle j is initialized with the circulation estimate at the end of analysis cycle j. Also, consider different forecast realizations resulting from perturbations in the forecast initial conditions of forecast cycle j that are consistent with the expected 4D-Var analysis error at the end of analysis cycle j. In addition, assume that the forecast surface forcing and open boundary conditions are also perturbed. The expected error covariance of the state vector relative to the unperturbed forecast is given by

 
formula

where (f)k is the tangent linear model at forecast time kΔt, linearized about the unperturbed forecast that is initialized using the unperturbed 4D-Var analysis, is the forecast perturbation derived from the analysis increment at the end of analysis cycle j, and , where is given by (10) and represents the state-vector forecast initial condition error covariance at the forecast start time, . Here, it is assumed that the forcing and boundary condition perturbations are drawn from distributions with the prior error covariances, and . In practice, however, these would be the error covariances of the forecast products used to drive the ocean forecasts. In general, would also include a component due to model errors, but in keeping with the strong constraint 4D-Var calculations of section 3, the model error covariance term is not included for simplicity. It is important to remember that the prior error covariance of the background for each data assimilation cycle is always , as prescribed in section 3, which in the present experiments is never updated using from (10). Therefore, the error estimates for each pair of consecutive analysis and forecast cycles can be considered to be independent.

Fig. 5.

(a) A schematic of adjacent analysis and forecast cycles. Shown are a 14-day forecast and a 7-day forecast initialized at the end of consecutive 4D-Var analysis cycles ending on days and , respectively, with both forecasts verifying on the same day, . The times , , , and are also referred to in the main text as k1Δt, k2Δt, k3Δt, and k4Δt, respectively. (b) As in (a), but illustrating different 14- and 7-day forecast realizations initialized from the analysis cycles ending on days and . The standard deviations (“spread”) of the 7- and 14-day forecast distributions are denoted as and , respectively.

Fig. 5.

(a) A schematic of adjacent analysis and forecast cycles. Shown are a 14-day forecast and a 7-day forecast initialized at the end of consecutive 4D-Var analysis cycles ending on days and , respectively, with both forecasts verifying on the same day, . The times , , , and are also referred to in the main text as k1Δt, k2Δt, k3Δt, and k4Δt, respectively. (b) As in (a), but illustrating different 14- and 7-day forecast realizations initialized from the analysis cycles ending on days and . The standard deviations (“spread”) of the 7- and 14-day forecast distributions are denoted as and , respectively.

Equation (13) is based on forecast realizations centered about a forecast initialized by the unperturbed analysis and is equivalent to the usual forecast ensemble employed at operational weather centers (e.g., Leith 1974; Molteni et al. 2006). Therefore, (13) also describes the expected forecast error covariance of a perfect model (i.e., no model error).

The covariance of any linear function of the forecast state vector xf(t) derived from all possible forecast realizations is given by , where k4Δt denotes the final forecast verification time (Fig. 5a). Figure 6 shows time series of for each index , which now are defined as the averages of the circulation over the last 48 h of 7-day forecasts initialized every week.

Fig. 6.

Time series (solid line) of for (a) , (b) , (c) , and (d) . In each case the contributions to of uncertainties in the initial conditions (line with open circles), surface forcing (line with closed circles), and boundary conditions (dashed line) are indicated, although in (b)–(d) the boundary condition contribution is negligible and is not shown. (e) Time series showing the contributions of uncertainties in the meridional wind stress (line with closed circles) and surface heat flux component (solid line) to in the case of .

Fig. 6.

Time series (solid line) of for (a) , (b) , (c) , and (d) . In each case the contributions to of uncertainties in the initial conditions (line with open circles), surface forcing (line with closed circles), and boundary conditions (dashed line) are indicated, although in (b)–(d) the boundary condition contribution is negligible and is not shown. (e) Time series showing the contributions of uncertainties in the meridional wind stress (line with closed circles) and surface heat flux component (solid line) to in the case of .

Also shown in Fig. 6 are time series of the contribution to of uncertainties in the initial conditions, surface forcing, and boundary conditions. With regard to the expected errors in the forecast transports, Fig. 6a reveals that the alongshore transport is controlled primarily by errors in the initial conditions, with errors in the surface forcing and open boundary conditions playing a minor role. The expected forecast errors for are generally largest during fall and winter. On the other hand, forecast error variance in the cross-shore transport is largely a result of uncertainties in the surface forcing (Fig. 6b). The biggest contributor by far in this case is errors in the alongshore wind stress (not shown), confirming the important role of Ekman transport. In this case uncertainties in the initial conditions account for about 25% of , and the forecast error variance is typically largest during the spring.

The expected forecast error variances for upper-ocean heat content, , and SST, (Figs. 6c and 6d), are also largest in the summer, particularly during the wind relaxation season in the case of (Fig. 6d, note the log scale), and are almost exclusively due to uncertainties in the surface forcing. The alongshore wind stress is generally the largest contributor to for during the summer, although uncertainties in the surface heat flux play a significant role also (Fig. 6e), especially during the winter. In the case of , errors in alongshore wind stress are the largest contributor to during summer, while at other times of the year initial condition errors are equally important. The large summertime contributions of surface forcing to forecast error are partly a reflection of the seasonal dependence imposed on the surface forcing prior error variances.

7. Predictability

To quantify the influence of data assimilation on the predictability of the circulation, consider again the overlapping analysis and forecast cycles j depicted in Fig. 5b. Consider now all of the possible realizations of 14- and 7-day forecasts initialized on days and , respectively, that verify on day . Each forecast realization is associated with uncertainties in the forecast initial conditions, surface forcing, and open boundary conditions drawn from a normal distribution with error covariance . The state vectors for a few of the 7- and 14-day forecast realizations are illustrated in Fig. 5b, and are denoted as and , respectively. The difference in the variance of the and realizations will be due solely to the assimilation of observations during the interval , assuming there are no model errors (or the effect of the model error is the same in both ensembles). If we denote the 7- and 14-day forecast variances of as and , respectively, then it can be shown that

 
formula

where and is the overlapping forecast period. The propagator 14 is linearized about the 14-day forecast , and a is the propagator for the jth analysis cycle spanning the interval . In the case of R4D-Var used here, a b 14, as discussed in section 2. Given that and are measures of the predictability of for and , then (14) quantifies the change in predictability due to assimilating observations. As noted in section 3, if 4D-Var is run to complete convergence, then , the true gain matrix, so using (3) then (14) reduces to as required according to (6).

Figure 7 shows a time series for each pair of consecutive analysis–forecast cycles of , which represents the percentage change in the forecast uncertainty due to assimilating observations. The case r = 100% represents the situation where the uncertainty in the 7-day forecast is reduced to zero by 4D-Var, while r = 0 corresponds to the case where 4D-Var yields no change in the uncertainty. When r < 0, this indicates that and that 4D-Var degrades the predictability of . Therefore, r > 0 indicates a positive impact of 4D-Var on the predictability of .

Fig. 7.

Time series of for (a) , (b) , (c) , and (d) . The colored bars show the contribution of each observation platform to r during each adjacent pair of analysis–forecast cycles. Key: SSH, Aviso SSH (red); SST, blended satellite SST (light green); Txbt, T from XBT casts (dark blue); Tctd, Sctd, T or S from CTD casts (brown and light blue); Targo, Sargo, T or S from Argo drifting floats (orange and yellow); and TOPP, T from tagged elephant seals (purple).

Fig. 7.

Time series of for (a) , (b) , (c) , and (d) . The colored bars show the contribution of each observation platform to r during each adjacent pair of analysis–forecast cycles. Key: SSH, Aviso SSH (red); SST, blended satellite SST (light green); Txbt, T from XBT casts (dark blue); Tctd, Sctd, T or S from CTD casts (brown and light blue); Targo, Sargo, T or S from Argo drifting floats (orange and yellow); and TOPP, T from tagged elephant seals (purple).

Figure 7 shows that r > 0 for all indicating that 4D-Var increases the predictability of . For , , and , the reduction in uncertainty is substantial, while for it is more modest. In the case of , the fraction r indicates a pronounced seasonal cycle in predictability, and data assimilation has the least impact on the uncertainty during wind relaxation periods, the same periods when forecast errors are typically largest (Fig. 6d). Since the vector in (14) resides in observation space, the fraction that each observation platform contributes to r can also be computed, and is indicated in Fig. 7 by the colored bars. For the alongshore transport, Fig. 7a shows that satellite SSH observations contribute most to the predictability, while for the cross-shore transport both SST and SSH contribute about equally on average. Not surprisingly, SST exerts the largest influence on the predictability of upper-ocean heat content and SST. However, for all , the in situ observations often exert considerable influence on r, despite being only about ~10% of the available data (M11b). We note here though that the impact of the observations on r will be influenced by the hypotheses embodied in , , and . Therefore, for example, overfitting of the observations or the presence of observations that, for any reason, are incompatible with the model may exert an overly large or undue influence on r.

As noted in the introduction, the use of the adjoint of a data assimilation system has recently found favor in NWP to quantify the impact of the observing system components on individual forecasts (e.g., Langland and Baker 2004; Gelaro and Zhu 2009). The focus of these studies has typically been on assessing the impact of each observation on a measure of forecast error, usually defined as the squared difference between the forecast and the verifying analysis at the same forecast lead time. In the approach used here, however, we consider the changes that arise in the spread of two theoretical, infinite-sized ensembles by using the adjoint of the 4D-Var algorithm to provide information about the expected errors that arise from perturbations in d drawn from the appropriate distribution. Thus, the information about the contribution of each observation to the predictability provided by our approach is quite different from observation impact information that has been considered in previous studies.

8. Consistency checks

In this section we present a series of consistency and reliability checks for the posterior and forecast errors derived in sections 57. It is important to reiterate that and given by (10) and (13) will be correct only if the prior error covariance and observation error covariance are correctly specified. The consistency checks presented here therefore provide an independent evaluation of the validity of the prior hypotheses about the errors in the background and the observations.

a. Prior minus posterior error

The expected covariance of the increments for any analysis cycle is given by . By the same token, for unbiased analyses , which shows that for any analysis cycle is drawn from a distribution with a variance . Therefore, if and are consistent with each other, the single realizations of for each R4D-Var cycle shown in Fig. 3 should lie within a distribution with a variance given by (11) and (12).

Figure 8 shows from each R4D-Var analysis cycle for the 7-day-average transport, heat content, and SST indices. Also shown in Fig. 8 are the standard deviation (std) ranges of ±1 std and ±2 std based on in Fig. 4. As noted in section 4, some of the s are biased, so Fig. 8 shows the bias-corrected values of (ab) for each index. For all indices, 35%–45% of the values of (ab) fall within ±1 std, and 69%–82% are within ±2 std, suggesting that for the most part, (ab) and are fairly consistent. For a large sample size, we might statistically expect ~95% of the values of (ab) to fall within ±2 std, so the 69%–82% range reported above indicates that some of the hypotheses that underpin and are incorrect and must be refined. Figure 8, however, does not offer any guidance on how to do this, only that there is some level of inconsistency.

Fig. 8.

Time series of (ab) (squares) from each R4D-Var analysis cycle for (a) 37N, (b) 500m, (c) HC, and (d) SST. Also shown in each case are time series of (solid lines) and (dashed lines).

Fig. 8.

Time series of (ab) (squares) from each R4D-Var analysis cycle for (a) 37N, (b) 500m, (c) HC, and (d) SST. Also shown in each case are time series of (solid lines) and (dashed lines).

b. Forecast error

The consistency and statistical reliability of the forecast error estimates can be assessed using ideas borrowed from probabilistic NWP (e.g., see Buizza et al. 2005). With this in mind, consider again the hypothetical infinite set of all possible forecast realizations in section 6 drawn from a normal distribution with error covariance f. If the forecasts are statistically reliable, then the true ocean state should belong to this same set of forecast realizations. In NWP the verifying analysis for the forecast time interval is usually taken as a surrogate for the truth and should, therefore, be a member of the same distribution as the forecasts that are initialized from the analysis at the end of the previous assimilation window. Therefore, on average, the distance between the mean forecast and the verifying analysis (e.g., the root-mean-square difference) should equal the average distance between the forecast mean and the forecast realizations (i.e., the forecast standard deviation). This situation is illustrated schematically in Fig. 9.

Fig. 9.

A schematic that illustrates an ensemble of 7-day forecasts xf of the state vector for the interval . Each ensemble member results from perturbing the forecast initial conditions, surface forcing, and boundary conditions with perturbations drawn from a distribution with covariance f. The unperturbed forecast is initialized at time using the posterior state vector estimate obtained at the end of analysis cycle j − 1 (see Fig. 5a). A few representative ensemble members are indicated, along with the unperturbed forecast, which is also equivalent to the ensemble mean in the case of the tangent linear assumption. Also indicated is the verifying analysis, which is the posterior state vector estimate obtained using 4D-Var during the interval , namely analysis cycle j (Fig. 5a). The difference between the unperturbed forecast (ensemble mean) and the verifying analysis is indicated as Δx, while the standard deviation of the ensemble about the ensemble mean (unperturbed forecast) is indicated as σf. For a statistically reliable forecast ensemble, . For each circulation index referred to in the main text, f is computed from the last 2 days of the unperturbed forecast (i.e., the interval ), while a is computed for the same interval using the verifying analysis.

Fig. 9.

A schematic that illustrates an ensemble of 7-day forecasts xf of the state vector for the interval . Each ensemble member results from perturbing the forecast initial conditions, surface forcing, and boundary conditions with perturbations drawn from a distribution with covariance f. The unperturbed forecast is initialized at time using the posterior state vector estimate obtained at the end of analysis cycle j − 1 (see Fig. 5a). A few representative ensemble members are indicated, along with the unperturbed forecast, which is also equivalent to the ensemble mean in the case of the tangent linear assumption. Also indicated is the verifying analysis, which is the posterior state vector estimate obtained using 4D-Var during the interval , namely analysis cycle j (Fig. 5a). The difference between the unperturbed forecast (ensemble mean) and the verifying analysis is indicated as Δx, while the standard deviation of the ensemble about the ensemble mean (unperturbed forecast) is indicated as σf. For a statistically reliable forecast ensemble, . For each circulation index referred to in the main text, f is computed from the last 2 days of the unperturbed forecast (i.e., the interval ), while a is computed for the same interval using the verifying analysis.

In section 6, the tangent linear assumption was employed to derive the forecast error covariance (13), in which case the mean forecast is simply the unperturbed forecast, xf (see Fig. 9). For each of the circulation indices, , we therefore expect the difference between the unperturbed forecast and the verifying analysis, (fa), to be drawn from a normal distribution with a standard deviation of based on (13). Since f also represents the mean forecast, it lies at the center of the forecast distribution, so significant departures of (fa) from the expected distribution indicates that the verifying analysis, a, is not a member of the forecast distribution, for which events the system is considered to be statistically unreliable.

Figure 10 shows the bias-corrected (fa) from each forecast cycle for the 2-day-average transport, heat content, and SST indices. Also shown in Fig. 10 are the ranges and of Fig. 6 for each cycle. The schematic in Fig. 9 indicates that for a statistically reliable forecast ensemble, we require , where Δx is the difference between the unperturbed forecast state vector and the state vector of the verifying analysis. For the linear scalar functions considered here, statistical reliability requires for each forecast cycle. The average of |(fa)| over all cycles is referred to as the mean absolute difference (MAD), while the average of over all cycles represents the mean standard deviation (MSTD), and both are indicated in Fig. 10 for each circulation index. For 37N and 500m, Fig. 10 indicates that 91% and 97% of the occurrences of (fa), respectively, fall within and , indicating that (fa) and are consistent most of the time. On the whole the system is statistically reliable for these indices since the verifying analysis a is statistically indistinguishable from a typical forecast realization, and MAD ~ MSTD. On the other hand, Fig. 10 reveals that only 40%–60% of the (fa) for heat content and SST fall within the range, indicating that (fa) and are inconsistent in these cases. For HC and SST, the forecast system is statistically unreliable because for most cycles the verifying analysis a is an outlier. Here, the unreliable nature of the forecast system is due to inappropriate choices of and . Given the large contribution of uncertainties in wind stress and heat flux to the forecast error in HC and SST (cf. Fig. 6), it seems likely that these components of are good candidates for improvement. However, model error has not been accounted for and may be the cause of much of the bias in HC and SST.

Fig. 10.

Time series of (fa) (squares) from each R4D-Var analysis cycle for (a) 37N, (b) 500m, (c) HC, and (d) SST. Also shown in each case are time series of (solid lines) and (dashed lines). In addition, the MAD given by over all cycles, and the MSTD given by over all cycles, are indicated. When MAD ≃ MSTD, the forecast system is considered to be statistically reliable.

Fig. 10.

Time series of (fa) (squares) from each R4D-Var analysis cycle for (a) 37N, (b) 500m, (c) HC, and (d) SST. Also shown in each case are time series of (solid lines) and (dashed lines). In addition, the MAD given by over all cycles, and the MSTD given by over all cycles, are indicated. When MAD ≃ MSTD, the forecast system is considered to be statistically reliable.

c. Posterior error

If sufficient observations are available, the analog of can be evaluated in observation space. If and are correctly prescribed, then following Desroziers et al. (2005)  , where ho is the equivalent of h in observation space. The scalar functions and are evaluated using the posterior estimate (xa) and the prior estimate (xb), respectively, evaluated at the observation points, and is computed from the observations (yo). The variance is the equivalent of (12) only now evaluated in observation space. Thus, provides a consistency check for computed from (12). Only a single realization of is available during each analysis cycle, but if it is consistent with (12), we would expect to be drawn from a distribution with the variance shown in Fig. 4.

In the present case, there are only sufficient observations available to estimate , the observation space equivalent of SST. Figure 11 shows for computed from each analysis cycle for the 7-day-average SST index. Due to the large range in , the absolute values are plotted so that a log scale can be used. Also shown in Fig. 11 are the variances and computed from (12) for each analysis cycle. Figure 11 reveals that the majority of the occurrences of lie outside the range of , indicating that the hypotheses embodied in and are in need of refinement.

Fig. 11.

Time series of (squares) from each R4D-Var analysis cycle for . Also shown are time series of (solid lines) and (dashed lines).

Fig. 11.

Time series of (squares) from each R4D-Var analysis cycle for . Also shown are time series of (solid lines) and (dashed lines).

d. Caveats

While the diagnostics presented here provide a useful indication of the consistency between the resulting analyses and forecast, and the assumptions made about the prior error covariance matrix and the observation error covariance matrix , they are probably not infallible. For example, the covariances and can be tuned to yield analysis errors and a cost function pattern of behavior that is closer to optimal, without the necessity of and being correct as demonstrated, for example, by Chapnik et al. (2004, 2006) and Desroziers et al. (2009).

9. Summary and conclusions

We have demonstrated how the adjoint of a 4D-Var data assimilation system can be used to compute estimates of the expected analysis and forecast error variances of linear functions of the ocean circulation. The use of adjoint 4D-Var in this way is a significant departure from methods currently used in meteorology and oceanography for estimating expected errors. Because the number of inner loops used during 4D-Var is typically much smaller than the number of observations, the posterior error covariance estimates based on a single, low-rank approximation of the true gain matrix do not effectively span the entire observation space, leading to an overestimate of the expected error. Alternatively, the adjoint of 4D-Var can be used to compute the sensitivity of the analysis increments to changes in the innovation vector resulting from changes in the control vector. By considering the expectation of all possible realizations of the analysis increments drawn from the appropriate normal distribution, the adjoint of 4D-Var can also be used to compute the posterior error covariance. Since the combination of the adjoint of 4D-Var and the expectation operator effectively spans the entire observation space, the resulting posterior error covariance estimates are more reliable than those based on . Due to the computational challenge involved, our approach is restricted at the present time to computing analysis (and forecast) error variance estimates for scalar functions of the state vector. However, the potential for using the adjoint of 4D-Var to estimate the full posterior and forecast error covariance matrices deserves further investigation. The expected error variance in quadratic scalar functions can also be estimated using the same approach if some additional assumptions are made. Quadratic functions will be discussed in a future publication.

In this study we have used the adjoint of ROMS 4D-Var to compute the analysis error variances of four indices that describe different aspects of the central California coast upwelling circulation pattern during a 2.5-yr sequence of 7-day 4D-Var data assimilation cycles. It was found that the expected error variances of the analyses from data assimilation are significantly lower than the prior error estimates. This is particularly true in the case of indices for which direct observations are available, as in the case of temperature, where the ratio of the posterior and the prior error variance reflects at least two orders of magnitude reduction.

The expected error variances of forecasts initialized from the 4D-Var analyses were also computed for the four circulation indices. A powerful aspect of the adjoint 4D-Var approach is that the contribution to the forecast error variance of uncertainties in the different components of the state-vector, surface forcing, and open boundary conditions can be readily computed. While some aspects of the forecast error variance behavior were expected, such as the large contribution of alongshore wind stress errors on the forecast skill of coastal SST, other findings were not anticipated, such as the apparent lack of impact of surface forcing uncertainties on the alongshore transport, or the negligible role played by errors in the open boundary conditions in all cases.

Another aspect of the adjoint 4D-Var approach with considerable utility is the ability to partition the predictability of the circulation indices (based on the standard deviation of all possible forecast realizations) into the contributions associated with each observation. This is in contrast to previous applications of data assimilation adjoints for computing the impact of observations on estimates of forecast error based on departures from a verifying analysis. In the system considered here, satellite data represent, by far, the largest fraction of available observations (M11b), and as such exert considerable influence on the forecast errors and predictability of the circulation. However, despite the relatively small fraction of subsurface observations, these data also sometimes exert a significant control on the forecast errors and predictability. The partitioning of the change in forecast uncertainty associated with each observation also allows the observations that impact the forecast skill the most to be identified, including those observations that are detrimental to the forecast.

A number of consistency checks were also presented, which provide an independent evaluation of the reliability of the analysis and forecast error estimates derived from the 4D-Var adjoint. On the whole, these checks indicate that the prior estimates of the background error covariance matrix and observation error covariance are in need of refinement. Nonetheless, despite the obvious and inevitable shortcomings of and , and the coarse resolution of the model used here, the utility of the adjoint 4D-Var approach is obvious, and we believe that the aforementioned conclusions are robust.

If a modular 4D-Var system is already in place, construction of the adjoint of 4D-Var is generally not difficult since all of the components required will already exist. In the case of ROMS, the inner-loop structure is composed of separate calls to the adjoint model, prior error covariance model, tangent linear model, and observation error covariance [see Fig. 3 in Moore et al. (2011a)], which facilitates calculation of the action of the stabilized representer matrix (T + ) on an input vector during R4D-Var. The minimization algorithm is also a separate module in ROMS, and the user may choose between the Lanczos formulation of the conjugate gradient algorithm or the minimum residual method (minres). Therefore, since the stabilized representer matrix is symmetric, the only additional model component required in ROMS for the calculations presented here was the adjoint of the minimization algorithm employed for identifying the minimum of the cost function, which is available for both the Lanczos algorithm and minres. The approach that we are advocating here is therefore practical for implementation in other models.

Acknowledgments

We are grateful for the continued support of the Office of Naval Research (N000140110209, N000140810556), and for support from the National Science Foundation (OCE-0628690) and the National Ocean Partnership Program (NA05NOS4731242). Any opinions, findings, and conclusions or recommendations expressed here are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors wish to thank Jim Doyle for the COAMPS surface forcing data, Dave Foley for the blended SST data, and Dan Costa and Patrick Robinson for the TOPP data that were all used in this study. Finally, we wish to thank the four anonymous reviewers, who provided valuable comments and feedback on earlier versions of this manuscript.

APPENDIX

Expected Analysis Error Covariance Estimates

In this  appendix we will compare the expected analysis error estimates based on the reduced-rank gain matrix and the adjoint dual 4D-Var operator (∂/∂d)T. The reduced-rank gain matrix is given by

 
formula

where m is the matrix of m orthonormal Lanczos vectors, qi, arising from the m inner loops of dual 4D-Var. Following Paige and Saunders (1975) and El Akkraoui and Gauthier (2010), the weight matrix , where m and m represent the LQ factorization of the m × m symmetric tridiagonal matrix m with diagonal elements and off-diagonal elements , where ai = qiδiqiγi−1qi−1 and is the stabilized representer matrix introduced in section 1. The matrix is an m × m diagonal matrix with leading diagonal elements equal to 1 except for the mth element, which is given by , where Lm,m is the mth diagonal element of m. Each Lanczos vector qi obeys the recurrence relation qi = γiqi+1 + δiqi + γi−1qi−1, where the first member of the sequence . Therefore, since m depends explicitly on d, the analysis increments can also be expressed as a function (d) as noted in section 2. Using the Lanczos vector expansion (A1), the tangent linearization of (d) can be expressed as

 
formula

where δm is the perturbation Lanczos vector matrix arising from the perturbations δd to the innovation vector, and .

The posterior error covariance matrix based on is given by . Equation (A1) shows that will at most map information from the m orthogonal directions of in observation space to the same m orthogonal directions of m via the matrix of weights m. Conversely the posterior covariance matrix given by (10) is found by considering all possible realizations of the innovation vector perturbations, δdl [drawn from a normal distribution with covariance (T + )], which in turn yield different realizations of the perturbation Lanczos vector matrices (δm)l that collectively span all directions in observation space. Based on (A2), (10) shows that will map the m orthogonal directions identified by the Lanczos vectors m into all of the directions in observation space (via the expectation operator) identified by all the realizations of the matrices (δm)l. Therefore, for the same choice of m, the operator (∂/∂d) provides a more complete mapping of a into observation space than . This is illustrated in Fig. A1, which shows the ratio of expected posterior error variance to the prior error variance, , versus m for each of the circulation indices introduced in section 4 computed using (∂/∂d) via (10) and via (A1) for the R4D-Var cycle spanning the period 5–12 June 2004. Figure A1 shows that for each , the ratio based on (∂/∂d) asymptotes to a constant value after about m = 60 inner loops, while based on is still decreasing even after m = 800 inner loops. This particular cycle was chosen for illustration because convergence of to its asymptotic value with increasing m is quite slow during this period. In general, however, as few as 25 inner loops were found to be adequate for reaching the asymptote of for the majority of the R4D-Var cycles (not shown).

Fig. A1.

The ratio vs the number of inner loops derived using ∂/∂d (circles) and using (squares) for (a) 37N, (b) 500m, (c) HC, and (d) SST. The horizontal dashed line denotes the value of reached using ∂/∂d after 60 inner loops.

Fig. A1.

The ratio vs the number of inner loops derived using ∂/∂d (circles) and using (squares) for (a) 37N, (b) 500m, (c) HC, and (d) SST. The horizontal dashed line denotes the value of reached using ∂/∂d after 60 inner loops.

REFERENCES

REFERENCES
Belo Pereira
,
M.
, and
L.
Berre
,
2006
:
The use of an ensemble approach to study the background error covariances in a global NWP model
.
Mon. Wea. Rev.
,
134
,
2466
2498
.
Bennett
,
A.
,
1992
:
Inverse Methods in Physical Oceanography. Cambridge University Press, 347 pp
.
Bennett
,
A.
,
2002
:
Inverse Modeling of the Ocean and Atmosphere. Cambridge University Press, 234 pp
.
Berre
,
L.
,
S.
Ştefnescu
, and
M. B.
Pereira
,
2006
:
The representation of the analysis effect in three error simulation techniques
.
Tellus
,
58A
,
196
209
.
Broquet
,
G.
,
C.
Edwards
,
A.
Moore
,
B.
Powell
,
M.
Veneziani
, and
J.
Doyle
,
2009a
:
Application of 4D-variational data assimilation to the California current system
.
Dyn. Atmos. Oceans
,
48
,
69
91
.
Broquet
,
G.
,
A.
Moore
,
H.
Arango
,
C.
Edwards
, and
B.
Powell
,
2009b
:
Ocean state and surface forcing correction using the ROMS-IS4DVAR data assimilation system. Mercator Ocean Quarterly Newsletter, No. 34, Ramonville St. Agne, France, 5–13
.
Broquet
,
G.
,
A.
Moore
,
H.
Arango
, and
C.
Edwards
,
2011
:
Corrections to ocean surface forcing in the California current system using 4D-variational data assimilation
.
Ocean Modell.
,
36
,
116
132
.
Buizza
,
R.
,
P.
Houtekamer
,
Z.
Toth
,
G.
Pellerin
,
M.
Wei
, and
Y.
Zhu
,
2005
:
A comparison of the ECMWF, MSC, and NCEP global ensemble prediction systems
.
Mon. Wea. Rev.
,
133
,
1076
1097
.
Chapman
,
D.
,
1985
:
Numerical treatment of cross-shelf open boundaries in a barotropic coastal ocean model
.
J. Phys. Oceanogr.
,
15
,
1060
1075
.
Chapnik
,
B.
,
G.
Desroziers
,
F.
Rabier
, and
O.
Talagrand
,
2004
:
Properties and first application of an error-statistics tuning method in variational assimilation
.
Quart. J. Roy. Meteor. Soc.
,
130
,
2253
2275
.
Chapnik
,
B.
,
G.
Desroziers
,
F.
Rabier
, and
O.
Talagrand
,
2006
:
Diagnosis and tuning of observational error in a quasi-operational data assimilation setting
.
Quart. J. Roy. Meteor. Soc.
,
132
,
543
565
.
Checkley
,
D.
, and
J.
Barth
,
2009
:
Patterns and process in the California current system
.
Prog. Oceanogr.
,
83
,
49
64
.
Chua
,
B.
,
L.
Xu
,
T.
Rosmond
, and
E.
Zaron
,
2009
:
Preconditioned representer-based variational data assimilation systems: Application to NAVDAS-AR. Data Assimilation for Atmospheric, Oceanic and Hydrologic Applications, S. K. Park and L. Xu, Eds., Springer-Verlag, 307–319
.
Cohn
,
S.
,
A. D.
Silva
,
J.
Guo
,
M.
Sienkiewicz
, and
D.
Lamich
,
1998
:
Assessing the effects of data selection with the DAO physical-space statistical analysis system
.
Mon. Wea. Rev.
,
126
,
2913
2926
.
Courtier
,
P.
,
J.-N.
Thépaut
, and
A.
Hollingsworth
,
1994
:
A strategy for operational implementation of 4D-Var using an incremental approach
.
Quart. J. Roy. Meteor. Soc.
,
120
,
1367
1388
.
Daget
,
N.
,
A.
Weaver
, and
M.
Balmaseda
,
2009
:
Ensemble estimation of background-error variances in a three-dimensional variational data assimilation system for the global ocean
.
Quart. J. Roy. Meteor. Soc.
,
135
,
1071
1094
.
Daley
,
R.
,
1991
:
Atmospheric Data Analysis. Cambridge University Press, 457 pp
.
Daley
,
R.
, and
E.
Barker
,
2001
:
NAVDAS: Formulation and diagnostics
.
Mon. Wea. Rev.
,
129
,
869
883
.
Desroziers
,
G.
,
L.
Berre
,
V.
Chabot
, and
P.
Poli
,
2005
:
Diagnosis of observation, background, and analysis-error statistics in observation space
.
Quart. J. Roy. Meteor. Soc.
,
131
,
3385
3396
.
Desroziers
,
G.
,
L.
Berre
,
V.
Chabot
, and
B.
Chapnik
,
2009
:
A posteriori diagnostics in an ensemble of perturbed analyses
.
Mon. Wea. Rev.
,
137
,
3420
3436
.
Doyle
,
J.
,
Q.
Jiang
,
Y.
Chao
, and
J.
Farrara
,
2009
:
High-resolution atmospheric modeling over the Monterey Bay during AOSN II
.
Deep-Sea Res. II
,
56
,
87
99
.
Egbert
,
G.
,
A.
Bennett
, and
M.
Foreman
,
1994
:
TOPEX/POSEIDON tides estimated using a global inverse method
.
J. Geophys. Res.
,
99
,
24 821
24 852
.
El Akkraoui
,
A.
, and
P.
Gauthier
,
2010
:
Convergence properties of the primal and dual forms of variational data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
136
,
107
115
.
Fairall
,
C.
,
E.
Bradley
,
J.
Godfrey
,
G.
Wick
,
J.
Ebson
, and
G.
Young
,
1996a
:
Cool-skin and warm layer effects on the sea surface temperature
.
J. Geophys. Res.
,
101
,
1295
1308
.
Fairall
,
C.
,
E.
Bradley
,
D.
Rogers
,
J.
Ebson
, and
G.
Young
,
1996b
:
Bulk parameterization of air–sea fluxes for Tropical Ocean Global Atmosphere Coupled Ocean–Atmosphere Response Experiment
.
J. Geophys. Res.
,
101
,
3747
3764
.
Fisher
,
M.
, and
P.
Courtier
,
1995
:
Estimating the covariance matrices of analysis and forecast error in variational data assimilation. ECMWF Tech. Memo. 220, Reading, United Kingdom, 28 pp
.
Flather
,
R.
,
1976
:
A tidal model of the northwest European continental shelf
.
Mem. Soc. Roy. Sci. Liege
,
6
(
10
),
141
164
.
Gelaro
,
R.
, and
Y.
Zhu
,
2009
:
Examination of observation impacts derives from observing system experiments (OSEs) and adjoint models
.
Tellus
,
61A
,
179
193
.
Gelaro
,
R.
,
T.
Rosmond
, and
R.
Daley
,
2002
:
Singular vector calculations with an analysis error variance metric
.
Mon. Wea. Rev.
,
130
,
1166
1186
.
Houtekamer
,
P. L.
,
L.
Lefaivre
,
J.
Derome
,
H.
Ritchie
, and
H. L.
Mitchell
,
1996
:
A system simulation approach to ensemble prediction
.
Mon. Wea. Rev.
,
124
,
1225
1242
.
Ide
,
K.
,
P.
Courtier
,
M.
Ghil
, and
A.
Lorenc
,
1997
:
Unified notation for data assimilation: Operational, sequential and variational
.
J. Meteor. Soc. Japan
,
75
,
181
189
.
Ingleby
,
B.
, and
M.
Huddleston
,
2007
:
Quality control of ocean temperature and salinity profiles—Historical and real-time data
.
J. Mar. Syst.
,
65
,
158
175
.
Isaksen
,
L.
,
M.
Bonavita
,
R.
Buizza
,
M.
Fisher
,
J.
Hasler
,
M.
Leutbecher
, and
L.
Raynaud
,
2010
:
Ensemble of data assimilations at ECMWF. ECMWF Tech. Memo. 636, Reading, United Kingdom, 48 pp
.
Lanczos
,
C.
,
1950
:
An iteration method for the solution of the eigenvalue problem of linear differential and integral operators
.
J. Res. Natl. Bur. Stand.
,
45
,
255
282
.
Langland
,
R.
, and
N.
Baker
,
2004
:
Estimation of observation impact using the NRL atmospheric variational data assimilation adjoint system
.
Tellus
,
56A
,
189
201
.
Leith
,
C.
,
1974
:
Theoretical skill of Monte Carlo forecasts
.
Mon. Wea. Rev.
,
102
,
409
418
.
Liu
,
W.
,
K.
Katsaros
, and
J.
Businger
,
1979
:
Bulk parameterization of the air–sea exchange of heat and water vapor including the molecular constraints at the interface
.
J. Atmos. Sci.
,
36
,
1722
1735
.
Molteni
,
F.
,
R.
Buizza
,
T. N.
Palmer
, and
T.
Petroliagis
,
2006
:
The ECMWF Ensemble Prediction System: Methodology and validation
.
Quart. J. Roy. Meteor. Soc.
,
122
,
73
119
.
Moore
,
A.
,
H.
Arango
,
G.
Broquet
,
B.
Powell
,
J.
Zavala-Garay
, and
A.
Weaver
,
2011a
:
The Regional Ocean Modeling system (ROMS) 4-dimensional variational data assimilation systems. I: System overview and formulation
.
Prog. Oceanogr.
,
91
,
34
49
.
Moore
,
A.
, and
Coauthors
,
2011b
:
The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems. II: Performance and application to the California current system
.
Prog. Oceanogr.
,
91
,
50
73
.
Paige
,
C.
, and
M.
Saunders
,
1975
:
Solution of sparse indefinite systems of linear equations
.
SIAM J. Numer. Anal.
,
12
,
617
629
.
Shchepetkin
,
A.
, and
J.
McWilliams
,
2005
:
The Regional Oceanic Modeling System (ROMS): A split explicit, free-surface, topography-following-coordinate oceanic model
.
Ocean Modell.
,
9
,
347
404
.
Tshimanga
,
J.
,
S.
Gratton
,
A.
Weaver
, and
A.
Sartenaer
,
2008
:
Limited-memory preconditioners with application to incremental variational data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
134
,
751
769
.
Weaver
,
A.
, and
P.
Courtier
,
2001
:
Correlation modelling on the sphere using a generalized diffusion equation
.
Quart. J. Roy. Meteor. Soc.
,
127
,
1815
1846
.
Weaver
,
A.
,
C.
Deltel
,
E.
Machu
,
S.
Ricci
, and
N.
Daget
,
2005
:
A multivariate balance operator for variational ocean data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
131
,
3605
3625
.
Wunsch
,
C.
, and
P.
Heimbach
,
2007
:
Practical global ocean state estimation
.
Physica D
,
230
,
197
208
.
Žagar
,
N.
,
E.
Andersson
, and
M.
Fisher
,
2005
:
Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors
.
Quart. J. Roy. Meteor. Soc.
,
131
,
987
1011
.
Zhu
,
Y.
, and
R.
Gelaro
,
2008
:
Observation sensitivity calculations using the adjoint of the Gridpoint Statistical Interpolation (GSI) analysis system
.
Mon. Wea. Rev.
,
136
,
335
351
.

Footnotes

1

The 4D-Var method proceeds by iteratively solving a sequence of linear, least squares minimization problems. The first member of the sequence is created via a series of inner loops linearized about the prior circulation. Subsequent members of the sequence are generated by updating the circulation about which the inner loops are linearized. This is done by rerunning the nonlinear circulation model using the newly computed increments from the last series of inner loops. The updates of the nonlinear model circulation in this way are referred to as outer loops.