## Jo Jo

If we introduce the shorthand notation s = (u,v,h) then Eqs. (14.1)-(14.3) can be written as st +As=0, where the subscript denotes differentiation with respect to time, and the operator A is given by:

The adjoint of (14.1)-(14.3) is defined by the Green's identity, which for the inner-product associated with the energy norm can be written as:

where s+= (u+,v+,h+) is a function in the space on which the adjoint A+ operates, and M= diag(H,H,g). To apply (14.6) to the shallow water equations, consider I = ¡0 ¡0 Hu+ x (14.1) + Hv+ x (14.2) + gh+ x (14.3)dxdy = 0. After integration by parts, it is easy to show that the adjoint of (14.1)—(14.3) is given by —s+ + A+s + = 0, where the negative time derivative indicates that time is reversed. The adjoint operator A+ is given by:

where the adjoint variables are periodic in x and satisfy the boundary conditions v+=0, at y = 0 and y = 1. In addition, s and s+ must satisfy the condition:

which is equivalent to time invariance of the inner-product {s+Ms}.

Comparing (14.5) and (14.7) shows that A+= -A with respect to the energy norm. In addition, the adjoint equation —s+ + A+s+ = 0 also satisfies the zero normal flow and periodic boundary conditions that are identical to those imposed on (14.1)— (14.3). Since A+=-A, the adjoint equation can also be written as — s+ — As + = 0.

Exercise2: Using I = ¡0 ¡0 Hu+ x (14.1) + Hv+ x (14.2) + gh+ x (14.3)dxdy = 0, derive the adjoint shallow water operator given by (14.7), and show that as a consequence of I = 0, the adjoint equation — s+ + A+s+ = 0 must satisfy (a) zero nor mal flow boundary conditions at y=0 and y = 1, and (b) the condition given by Eq. (14.8).

Wave solutions of st +As=0 given by (14.1)-(14.3) take the form of eastward and westward propagating inertia-gravity waves and Rossby waves (Gill 1982). Similarly wave solutions exist for the adjoint equation -s+ - As + = 0 (recall A+= -A) but with opposite phase and group velocities because of the reversal of time. For example, long wavelength Rossby waves that carry energy westward in the shallow water equations s t +As=0 will carry energy eastward in the adjoint equations -s+ - As + = 0.

The adjoint equation can also be written as s+ + As+ = 0 showing that it is mathematically equivalent to Eqs. (14.1)-(14.3). Therefore, for the same initial conditions the solutions s=(u,v,h) and s+=(u+,v+,h+) will be identical, in which case Eq. (14.8) is an expression of energy conservation.

14.2.3.2 The Linear Shallow Water Equations in the Presence of a Mean Circulation

Consider now the case of linear waves in the same periodic rectangular ocean which is now in motion with a mean geostrophic circulation u = ( - g / f ) dh/dy. The linearized shallow water equations in this case are:

where H = H + h, and h is the sea surface displacement associated with the geo-strophic circulation.

Using the same compact form as before, we can express (14.9)-(14.11) as s t+As=0, where the operator A is now given by:

Applying the Green's identity and using the energy norm, the adjoint equation is of the form -s+ + A+s+ = 0 where now:

As before, Eq. (14.8) is a necessary condition, and the adjoint variables are periodic in x and satisfy the boundary condition v+ = 0, at y=0 and y = 1. In this case, the y-gradient of the mean circulation u can act as a source of energy for linear waves, and is the familiar source of barotropic instability if the waves can undergo sustained exponential growth. In this case, solutions of st+As = 0 and s+ - A+s + = 0 with the same initial conditions will no longer be identical, and energy is no longer conserved. Sustained exponential growth will occur if the potential vorticity gradient df/dy — d2u/d2y changes sign anywhere within the channel (Pedlosky 1987).

Exercise 3: Using I = f01 f01 Hu+ x (14.9) + Hv+ x (14.10) + gh+ x (14.11)dxdy = 0, derive the adjoint shallow water operator given by (14.13), and show that as a consequence of I = 0, the adjoint equation -s+ + A+s + = 0 must satisfy (a) zero normal flow boundary conditions y=0 and y = 1, and (b) the condition given by Eq. (14.8).

### 14.3 Variational Data Assimilation 14.3.1 Notation

Before proceeding to describe the important role of adjoint operators in variational data assimilation methods, it is necessary to introduce some notation. Ocean models solve the discrete equations of motion on grids in both space and time. There are a wide range of ocean models available (see chapters by Barnier and Chassignet), and many use different kinds of grid configurations and coordinate systems (e.g., staggered grids, terrain following vertical coordinates, isopyc-nal vertical coordinates, etc). Nonetheless, all models solve for a standard set of prognostic variables, typically temperature, salinity, velocity, and free surface displacement, and can be expressed in a generic symbolic form using the concept of a state-vector. To this end we introduce a state-vector x(t.) which represents a vector of all grid point values of the prognostic state variables in space at a given time t.. The ocean state x(t.) will depend on the state at some earlier time x(ti-1) and upon the ocean surface forcing f(ti) and boundary conditions b(ti) over the time interval [ti-1,ti]. The discrete equations of motion will in general be nonlinear, and can be represented by the discrete nonlinear operator M. Thus the time evolution of the ocean state by an ocean model can be expressed in a convenient and compact form as:

where M(ti,ti-1) denotes a forward integration of the nonlinear ocean model from time ti-1 to ti and for convenience f(t) and b(ti) denote the forcing and boundary conditions over the entire interval [ti-1,ti]. This notation is fairly standard in numerical weather prediction and ocean modeling (Ide et al. 1997).

### 14.3.2 The Incremental Formulation

As discussed in the chapters by Zaron and Brasseur, the aim of data assimilation is to construct an estimate of the ocean circulation by combining prior information from an ocean model with observations. The solution of an ocean model over the interval t = [t0,tN] is uniquely determined by the initial conditions, x(t0), the surface forcing, f(t), and the boundary conditions, b(t), collectively referred to as control variables. Therefore, prior estimates of all control variables are required, and will be denoted xb(t0), fb(t) and bb(t) respectively. Our hypothesis about the uncertainty associated with each prior is embodied in the prior error covariance matrices for the initial conditions, B, the surface forcing, Bf, and boundary conditions, Bb. For convenience, we will denote by D the block diagonal covariance matrix comprised of B, Bf, and Bb, namely D = diag(B, Bf, Bb). Similarly, we will denote by y the vector comprised of all observations in the interval t = [t0,tN] with the associated observation error covariance matrix R. The goal of data assimilation is then to identify z = (xT(io), fT(io), fT(ii),..., fT(iw), bT(io), bT(ii),..., bT(tN ))T, the so-called control vector, that maximizes the conditional probability p(z|y) a e-jNL, where:

and 9 is the vector of the model equivalent of the observations at the observation times and locations. Each element of 9 is of the form H (x(t)) where H is the operator that transforms or interpolates the state-vector x(t) to the observation points at time t The scalar JNL is called the penalty function or cost function.

Since the model M and observation operator H are in general nonlinear, JNL is a non-quadratic form so it may be convex and there may be multiple values of z that minimize p(z|y). Therefore, it is common to linearize M and H by considering small increments bz to the prior zb, so that z=zb + <5z (Courtier et al. 1994), where Zb = (xT(io), fbT(io), fb (ii), • • •,bT(io),bT(ii),... )T is the vector of priors. The assumption underlying this approximation is that zb does not lie too far from the true state of the ocean, in which case the goal of data assimilation then becomes one of finding 8z = (SxT(t0), SfT(t0), SfT(t1),..., SbT(t0), SbT(t1),... )T, the vector of control variable increments, that minimizes the linearized form of (14.15), namely:

J(Sz) = 2SzTD-1Sz + 1 (GSz - d)TR-1 (GSz - d) (14.16)

where the matrix G is the operator that maps the control variable increments to the observation points. The vector d=y - 9b is called the innovation vector, where 9b is the vector (H(xb(/.))). The increment bza that minimizes J corresponds to the maximum value of p(5z|y), and is often called a maximum likelihood estimate. The maximum likelihood ocean state estimate is then given by the analysis or so-called posterior za=zb + bza. In the event that the prior hypotheses embodied in D and R are correct, then the theoretical minimum value of the cost/penalty function is J = N , /2, half the total number of observations (Bennett 2002).

The prior circulation estimate xb(t) during the interval t = [t0,tN] is assumed to be a solution of the model Eq. (14.14) forced by the prior fb(t) and subject to the prior boundary conditions bb(t). Under the assumption that the prior is already a good estimate of the circulation, the increments dz will be small compared to zb, in which case a good approximation for dx(t) will be the first-order Taylor expansion of (14.14), namely:

where M(ti,ti-1) represents the linearization of M in (14.14) about the time evolving prior xb(t), and Su(ti-1) = (SxT(ti-1), SiT(ti), SbT(ti))T. Equation (14.17) is referred to as the tangent linear model, since solutions dx(/.) are locally tangent to the solution xb(t.) of (14.14). The operator G in (14.16) is a convolution in time of the tangent linear model M and H (the linearization of the observation operator H), and represents solutions of (14.17) evaluated or mapped to the observation points.

The increment dza corresponding to the most likely state estimate satisfies the condition dJfdSz = 0, and is given by dza=Kd where K is called the gain matrix and given by:

Equivalently, the gain matrix can also be written as:

In both (14.18) and (14.19), evaluation of the gain matrix involves a matrix inverse. In (14.18), the matrix to be inverted is (D-1 + GTR-1 G) and has the dimension of dz which will generally be greater than the number of model grid points and may be very large and a challenge to invert. The dimension of dz is Nm=(Nx+Nf+Nb) where Nx, is the dimension of x, and Nf and Nb are the dimensions of f, and b, respectively, multiplied by the number of model time steps N in the interval t = [t0,tN]. In addition, the expression in (14.18) in parentheses involves D-1 which may also be difficult to evaluate. Eq. (14.18) is often referred to as the primal form of the gain matrix.

Alternatively the matrix to be inverted in (14.19) is (GDGT + R) which has a dimension equal to that of y, the number of observations, Nobs. In general, Nobs ^ Nm, so (14.19) may be more convenient to use than (14.18). Equation (14.19) is often referred to as the dual form of the gain matrix. In practice, both (14.18) and (14.19) are used for ocean data assimilation and in either case the most likely state of the ocean circulation is given by za=zb+Kd.

Exercise 4: Prove that za=zb+Kd when dJfdSz = 0, and that K is given by (14.18).

Exercise 5: Using the identity (A+B)-1=A-1(A-1 + B-1)-1B-1, show that K can also be expressed in the dual form (14.19).

Regardless of whether the primal or dual form of K is used, the matrix inverse in (14.18) and (14.19) is never explicitly evaluated. Instead, an equivalent system of linear equations is solved, and the adjoint GT of the tangent linear model sampled at the observation points, G, plays a crucial role in this process. The adjoint of the tangent linear model can be expressed symbolically as:

where the reversed order of the time arguments of MT compared to M in (14.17) indicates that integration is backwards in time, as in the shallow water examples of Sect. 14.2.3. The matrix GT is a time convolution of the adjoint model MT and HT, the adjoint of the linearized observation operator. Since identification of the most likely increment Kd is equivalent to identifying the condition dJ/dSz = 0, the resulting data assimilation methods are collectively referred to as 4-dimensional variational (4D-Var) data assimilation, where the four dimensions are space and time.

As in Sect. 14.2.1, there are two fundamental spaces, the primal space of dimension Nm, and dual space of dimension Nob. The tangent linear operator G maps a vector from primal space to dual space, while the adjoint of the tangent linear operator GT maps a vector from dual space to primal space.

### 14.3.3 Primal Space 4D-Var

It is important to realize that while (14.18) and (14.19) are written in terms of matrix products, the matrices involved are never explicitly computed, and all matrix manipulations are performed using models, including models for D and R. This leads to very useful iterative algorithms that can be used to identify the minimum of J(5z) or JNL(z). With this in mind, consider the derivative of J(5z) in (14.16) with respect to ¿z:

Equation (14.21) shows that the gradient of the cost/penalty function with respect to ¿z can by evaluated by (1) running the tangent linear model subject to ¿z to evaluate G£z, the tangent linear model solution sampled at the observation points, (2) running the adjoint of the tangent linear model, GT, forced by R-1 ^¿z-d), the weighted difference between G¿z and the innovation vector, d, and (3) adding D-^z to the result of (2). Therefore, the cost/penalty function (14.16) can be minimized itera-tively as follows:

1. Choose an initial starting value of ¿z.

2. Run the nonlinear model (14.14) with the prior initial conditions, xb(t0), prior forcing, fb(t), and prior boundary conditions, bb(t), and compute the prior circulation estimate xb(t) for the interval t = [t0,tN].

3. Run the tangent linear model (14.17) linearized about xb(t) from step 2 for the interval t = [t0,tN], compute J^z) from (14.16), and compute and save R-1(G¿z-d).

4. Run the adjoint of the tangent linear model (14.20) backwards in time linearized about xb( t) from step 2 and forced by R-1(GSz-d) from step 3 for the interval t = [twto].

5. Add to the adjoint solution from step 4 the vector D-1Sz to yield dJ/dSz according to (14.21).

6. Using the cost function gradient dJ/dSz from step 5, use a conjugate gradient method to identify a new dz that will reduce the value of J resulting from a subsequent run of the tangent linear model as in step 3.

7. Using the new Sz from step 6, repeat steps 3-6 until the minimum of J has been identified.

At the minimum of J, the gradient dJ/dSz = 0 and the iterative procedure described by steps 2-7 is equivalent to evaluating Kd using the primal form of the gain matrix in (14.18).

### 14.3.4 Dual Space 4D-Var

In the dual formulation of the optimal increment Kd, the matrix in parentheses in (14.19) is inverted by introducing an intermediate variable w = (GDGT + R)-1d, where Kd=DGTw. In practice, w is identified by solving the linear system (GDGt + R)w = d using an iterative conjugate gradient method to minimize the function I = l/2wT(GDGT + R)w - wT d, and w plays the role of a generating function as discussed in Sect. 14.2.1. The increment Kd is then identified according to DGTw. The following iteration algorithm is typical of those commonly used:

### 1. Choose an initial starting value of w.

2. Run the nonlinear model (14.14) with the prior initial conditions, xb(t0), prior forcing, fb(t), and prior boundary conditions, bb(t), and compute the prior circulation estimate xb(t) for the interval t = [t0,tN].

3. Run the adjoint of the tangent linear model (14.20) backwards in time linearized about xb(t) from step 2 and forced by w for the interval t = [tN,t0] to yield GTw.

4. Apply the prior covariance D to the adjoint model solution from step 3 at t = 0 to yield DGTw.

5. Run the tangent linear model (14.17) linearized about xb(t) from step 2 for the interval t = [t0,tN] using the result of step 4 as the initial condition to yield GDGTw.

6. Add Rw to the result of step 5, and evaluate the gradient d 1/dw = (GDGt + R)w - d.

7. Using the gradient dI/dw from step 6, use a conjugate gradient method to identify a new w that will reduce the value of I resulting from a subsequent reevaluation of steps 3-5.

8. Using the new w from step 7, repeat steps 3-7 until the minimum of I has been identified.

9. Having identified the w=wa that minimizes I, evaluate the increment Sz = Kd=DGTw by repeating steps 3 and 4 using w .

Fig. 14.1 A schematic illustrating the action of the adjoint, prior error covariance, and tangent linear operators on a 8 function located at the site of a single observation in a steady zonal shear flow represented by the blue arrows. a The initial location of the 8 function in dual space. b The action of the adjoint operator GT which propagates the 8 function backwards in time "upstream" and maps it into primal space. c The prior error covariance matrix smoothes GT8 in primal space. d The final action of tangent linear operator G propagates the smoothed field from step (c) forward in time and maps the field back to dual space, indicated by the red open circle located at the observation point u a b c

The connection between primal and dual space, and the role that G and GT play in transforming between one space and the other is probably best illustrated by considering the operations represented by steps 3, 4 and 5 when applied to a single observation. Consider a ¿-function at the time and location of the single observation, in which case steps 3-5 yield GDG1^. The sequence of operations that lead to this result are illustrated schematically in Fig. 14.1 for the case of a mean geostrophic flow in the form of a jet in the zonal channel considered in Sect. 14.2.3.2. In this example, advection is the dominant dynamical process, although there will be some wave propagation as well as illustrated in Fig. 14.1b. Hence the actions of GT and G primarily advect information upstream and downstream respectively.

### 14.3.5 Computation of za

Having identified the increment control vector <za using either the primal or dual form of 4D-Var, it then remains to compute the most likely circulation estimate xa(t)=xb(t) + <5xa(t) over the interval t = [t0,tN]. Two approaches are generally used: (1) using the non-linear model (14.14) to advance the circulation x(t0) = xb(t0) + <xa(t0) forward in time using fb(t)+<fa(t) and bb(t)+<ba(t), or (2) using the tangent linear model (14.17) forced by <f (t) and subject to <b (t) to advance the increments bx(t) in time. Approach (1) is generally used in primal space applications of 4D-Var (Courtier et al. 1994), while both (1) and (2) are used in dual space formulations (Da Silva et al. 1995; Egbert et al. 1994). To the author's knowledge there are no primal formulations that use (2). Of course in the case of a linear model, (1) and (2) are equivalent.

In the dual formulation, each element of the innovation vector d is assumed to be a linear combination of the elements of the state-vector increment bx, according to the tangent linear operator G, and it is the generating function w of Sect. 14.3.4 that identifies the activated part of primal space into which d maps. All linear functions of bx are known as the dual of bx, and possess the important property that the tangent linear model equivalent of each element d. of d can be expressed as rj Sx, the so-called Riesz representation theorem. The time dependent vectors r. are called representer functions, and approach (2) in the dual formulation is equivalent to expressing the increment bx (t) as Sc, where S is a N x N , matrix, and each column it o a\ / ? m obs '

of S is a representer function r. The vector c is comprised of the weights assigned to each of the Nobs representers (Bennett 2002). In the schematic of Fig. 14.1d, the pattern of colored contours is a representer function, and the open red circle is GDGTb, the representer function sampled at the observation point.

### 14.3.6 Strong Constraint Versus Weak Constraint

In the formulations of 4D-Var presented so far, it has been implicitly assumed that the most likely circulation estimate xa(t) is an exact solution of the nonlinear model Eq. (14.14). This is tantamount to assuming that the model is perfect and free of errors, and the increments bza(t) that minimize the cost/penalty function Jin (14.16) are said to be subject to the "strong constraint" imposed by model dynamics (Sasaki 1970). Of course, all models possess errors and uncertainties, and to account for these it is necessary to augment the control vector of increments bz so that:

5z = (SxT(to),sfT(to),sfT(to,..., sbT(to), 5bT(ti),..., n(to), n(fi),.. .)T (14.22)

where n(t) represents the corrections for model error at each grid point and time step. The prior for model error is assumed to be 0 with an associated error covari-ance matrix Q. In the presence of model error, the development of 4D-Var in primal and dual space of Sects. 14.3.1-14.3.5 is unchanged, except that now the prior error covariance matrix is given by the block diagonal matrix D=diag(B,Bf,Bb,Q). The most likely circulation estimate xa(t) in this case is no longer an exact solution of the non-linear model Eq. (14.14), and the increments bza(t) that minimize the cost/ penalty function J(bz) in (14.16) are said to be subject to the "weak constraint" imposed by model dynamics (Sasaki 1970).

Under the weak constraint, the dimension, N , of primal space increases by N N, and in general the weak constraint 4D-Var problem becomes intractable. However, the dimension of dual space is unchanged by the imposition of the weak constraint, so weak constraint 4D-Var is often performed using the dual formulation.

### 14.3.7 Inner- and Outer-Loops

Following the incremental formulation of 4D-Var described in Sect. 14.3.2, the most likely circulation is that which minimizes the cost function J(ôz) given by (14.16). However, the assumption underlying (14.16) is that the prior is close to the true circulation. This of course is a big assumption, and in all likelihood will not always be true, if ever. Therefore, it is preferable to identify the most likely circulation estimate that minimizes instead the cost function J^z) in (14.15). In practice, JNL(z) is a non-quadratic function of z because the state-vector is a solution of the nonlinear model (14.14). As a result, JNL(z) may possess multiple minima, and a global minimum corresponding to the most likely circulation may be difficult to identify. However, a common technique for identifying the minima of (14.15) is to solve a sequence of linear minimizations of the form (14.16), where each member of the sequence is referred to as an "outer-loop." During each outer-loop, J(<z) given by (14.16) is minimized using the iterative algorithms of Sect. 14.3.3 or Sect. 14.3.4, and each iteration is called an "inner-loop."

During the first outer-loop, the tangent linear model (14.17) and adjoint model (14.20) are linearized about the prior circulation estimate xb( t) over the interval t = [t0,tN]. At the end of the first outer-loop, the circulation estimate is updated using approach (1) or (2) of Sect. 14.3.5, to yield the state vector, Xj(t)=xb(t) + <5xj(t), forcing, fj(t)=fb(t) + <5fj(t), boundary conditions, bj(t) = bb(t) + <5bj(t), and in the weak constraint case the corrections for model error, ^(t), where the subscript "1" refers to the first outer-loop. During the second outer-loop, the tangent linear and adjoint models are linearized about Xj(t). Repeating the sequence of outer-loops, it is easy to see that during outer-loop n, the tangent linear and adjoint models are linearized about xn-1(t), where xn-1(t) is the updated circulation estimate forced by fn-1(t), and subject to bn-1(t) and nn-1(t). It is important to note, however, that the innovation vector never changes, and d=y - (H.(xb(t.))) is always computed using the prior circulation estimate xb(t.).

The method used to compute the circulation estimate xn(t) at the end of each outer-loop varies. The most common approach is to use the non-linear model as in approach (1) of Sect. 14.3.5 to advance Xn(to) = xb(to) + 51"= ^(to) forward in time. However, in the dual formulation in which the method of representers is used, approach (2) is modified, and the finite-amplitude tangent linear model is used to advance xn(t0) forward in time. The finite-amplitude tangent linear model can be expressed symbolically as:

Xn(ti ) = M (ti, ti-i)(x„-i(ti-i), fn-i(ti ), b„-i(ii )) + M„_i(ti, ti-i)(gn(ti-l) - gn-l(ti-l))

where gn(ti) = (x'T(ti), fTfc), bTn(ti ))T, and Mb-1 denotes the tangent linear model linearized about xB-1( t). During the first outer-loop, the solution of (14.23) reduces to the sum of the prior, xb(t), and the tangent linear model solution (14.17). The representer approach to 4D-Var in dual space using the algorithm of Sect. 14.3.4 for the inner-loops, and updating the circulation estimates in the outer-loops us-

ing (14.23) is equivalent to minimizing JNL(z) in (14.15) by solving the non-linear Euler-Lagrange equations. Full details of this approach are beyond the scope of this article but can be found in Bennett (2002).

### 14.4 Examples of 4D-Var for the California Current

We will present here some illustrative examples of 4D-Var using both the primal and dual formulations applied to the California Current circulation using the Regional Ocean Modeling System (ROMS).

14.4.1 The Regional Ocean Modeling System (ROMS)

ROMS is a primitive equation ocean model that has gained considerable popularity in recent years because of the great flexibility that it affords for modeling different regions of the world ocean. ROMS uses a curvilinear orthogonal coordinate system in the horizontal, and terrain-following coordinates in the vertical, both of which allow for increased resolution in regions where it is most needed (e.g., in regions of complex topography and bathymetry, in shallow water, and near the ocean surface). ROMS is a hydrostatic model, and employs a wide range of user-controlled options for the numerics and physical parameterizations, as well as a range of options for prescribing open boundary conditions. A detailed description of ROMS is beyond the scope of this article, and the reader is referred to Shchepetkin and McWilliams (2005) and Haidvogel et al. (2000) for more information. ROMS is a community ocean model and is freely available from http://www.myroms.org.

### 14.4.2 ROMS 4D-Var

There are many practical aspects of 4D-Var that will not be discussed here, but which nonetheless are very important. The main features of the ROMS 4D-Var system are listed below, along with references that provide more information. In addition, the ROMS 4D-Var system is described in detail by Moore et al. (2011a, b, c). The main features and attributes of ROMS 4D-Var can be summarized as follows:

1. Incremental formulation (Courtier et al. 1994).

2. Primal and dual formulations (Courtier 1997).

3. Primal formulation is referred to as Incremental 4D-Var, hereafter I4D-Var.

4. Dual formulations following the Physical-space Statistical Analysis System of Da Silva et al. (1995), hereafter 4D-PSAS, and the indirect representer method of Egbert et al. (1994), hereafter R4D-Var.

5. Strong constraint formulation for both primal and dual formulations.

### 6. Weak constraint for dual formulations only.

7. Inner-loops use a preconditioned Lanczos formulation of the conjugate method (Golub and van Loan 1989; Lorenc 2003; Fisher and Courtier 1995; Tshimanga et al. 2008).

8. Prior covariances D are modeled using a pseudo-heat diffusion equation approach (Derber and Bouttier 1999; Weaver and Courtier 2001), and a multivariate balance operator (Weaver et al. 2005).

9. MPI parallel architecture.

### 14.4.3 The California Current System

The California Current System (CCS) is a prototype eastern boundary current regime that is dominated by mesoscale eddies, and characterized by a pronounced seasonal cycle of coastal upwelling and primary productivity. A comprehensive review of the CCS can be found in Hickey (1998).

ROMS has been configured for the CCS (hereafter referred to as ROMS-CCS) and spans the domain 116W-134W, 30N-48N shown in Fig. 14.2. Several configu-

Fig. 14.2 The ROMS-CCS domain and bathymetry in meters on a grid with 10 km horizontal resolution 5000 4500 4000 3500 3000 2500 2000 1500 1000 500

Fig. 14.2 The ROMS-CCS domain and bathymetry in meters on a grid with 10 km horizontal resolution

### 5000 4500 4000 3500 3000 2500 2000 1500 1000 500

rations of ROMS-CCS exist, with horizontal resolutions ranging from 3 to 30 km, and with 30-42 levels in the vertical. ROMS-CCS is forced with near surface air data from the Coupled Ocean Atmosphere Mesoscale Prediction System (COAMPS) of Doyle et al. (2009) which are converted to surface wind stress and surface fluxes of heat and freshwater using the bulk flux formulations of Fairall et al. (1996a, b) and Liu et al. (1979). The model domain has open boundaries at the northern, western, and southern edges, and at these boundaries the temperature, salinity, and velocity fields are prescribed using data from the Estimating the Circulation and Climate of the Ocean project (ECCO) (Wunsch and Heimbach 2007) which is a global ocean data assimilation product at 1° resolution. The radiation conditions of Chapman (1985) and Flather (1976) are also imposed at the open boundaries on the free surface and vertically integrated velocity respectively. ROMS-CCS produces very good simulations of the CCS circulation as documented by Veneziani et al. (2009).

### 14.4.4 ROMS-CCS 4D-Var Configuration

ROMS-CCS has been used extensively for 4D-Var as described by Broquet et al. (2009a, b, 2011) and Moore et al. (2011b, c). Some example 4D-Var calculations using ROMS-CCS are presented here to illustrate the ideas and concepts introduced in Sects. 14.2 and 14.3.

Recall that for data assimilation we require prior estimates for the elements of the control vector z, namely for the initial conditions, xb(t0), the surface forcing, fb(t), and the boundary conditions, bb(t), for the interval t = [t0,tN]. In the case of weak constraint 4D-Var, the prior for the corrections for model error is the null vector 0(t). ROMS-CCS 4D-Var is usually run sequentially using data assimilation windows in time that span time intervals tN-t0 that are typically 4-14 days long. Each data assimilation window is referred to as a "cycle," and the initial condition prior, xb(t0), is the best circulation estimate from the end of the previous cycle. The priors for surface forcing, fb(t), are the surface fluxes derived from the COAMPS air data, and the priors for the boundary conditions, bb(t), are the open boundary data from ECCO. b

In addition to the prior fields, prior error covariance matrices are also required, namely, B for the initial conditions, Bf for the surface forcing, Bb for the boundary conditions, and Q for the model errors in the case of weak constraint 4D-Var. In general, the prior error covariance matrices vary in time, but in practice they are assumed to be time invariant. Following Weaver et al. (2005), each of the prior error covariance matrices are factorized into a block diagonal, univariate correlation matrix, a diagonal matrix of standard deviations, and a multivariate dynamical balance operator. The univariate correlation matrix is modeled as the solution of a pseudo-heat diffusion equation (Weaver and Courtier 2001). This is an involved procedure, and the prescription of each prior error covariance matrix is beyond the scope of this presentation. However, a full description for ROMS-CCS can be found in Broquet et al. (2009). For the present purpose, it suffices to say that the prior er ror covariance matrices can be prescribed for each of the prior components of the control vector zb.

Observations from various platforms are assimilated into ROMS-CCS, including satellite derived sea surface temperature (SST) and sea surface height (SSH), and sub-surface hydrographic measurements of temperature, T, and salinity, S, collected from shipboard CTDs, XBTs, and drifting Argo profiling floats. The data used were extracted from the rigorously quality controlled global ocean data archive of Ingleby and Huddleston (2007). The observation error covariance matrix, R, is assumed to be time invariant and diagonal (i.e. spatially and temporally uncorrelated observation errors), and the error variances are instrument and variable dependent. The elements of R are a reflection of several sources of error: instrument error, interpolation errors introduced by the observation operator H of Sect. 14.3.2, and errors of representativeness. The largest source of error is typically the error of representativeness, which is a measure of the uncertainty associated with the ability of a single ocean observation to describe the circulation in a single ocean model grid cell. The following observation error standard deviations were used in ROMS-CCS: ~2 cm for SSH; ~0.4C for SST; ~0.1C for T; ~0.01 for S.

### 14.4.4.1 Primal Versus Dual 4D-Var

The first example considered will contrast the performance of the primal and dual formulations of ROMS 4D-Var in the CCS, and results are shown from the ROMS-CCS configuration with a horizontal resolution of 30 km and 30 levels in the vertical. The cost function from two strong constraint calculations using 1 outer-loop and 75 inner-loops, is shown Fig. 14.3 for a 4 day assimilation window spanning the period 1-4 July, 2000, during which time there were ~1 x 104 observations available, most in the form of satellite data.

In one calculation the primal formulation I4D-Var was used, while in the second calculation the dual formulation R4D-Var was employed. Figure 14.3 shows that in both cases the cost function decreases with an increasing number of inner-loop iterations. In the case of I4D-Var, J(bz) decreases monotonically, and asymptotes to a near constant value after about 40 or so inner-loops. R4D-Var on the other hand exhibits very different behavior in which J(bz) undergoes quite large fluctuations until after about 60 inner-loops it asymptotes to a similar near constant value as that reached by I4D-Var. The difference in behavior of J(bz) in the primal and dual formulations of 4D-Var is associated with the different approach employed to identify the cost function minimum. In the primal case, J(bz) is minimized directly as described in Sect. 14.3.3, while in the dual formulation J(bz) is minimized indirectly using a generating function w according to Sect. 14.3.4. El Akkraoui and Gauthier (2009) have explored the different behavior of the primal and dual formulations and suggest some effective remedies for increasing the rate of convergence of J(bz) in the dual case.

Nonetheless, despite the remarkable difference between the primal and dual formulations in the rate of convergence of J( bz) to its asymptotic minimum value, Number of inner-loops

Fig. 14.3 Log10(J(5z)) (from (14.16)) versus the number of inner-loops for three 4D-Var assimilation calculations for the period 1-4 July, 2000 for the case of 1 outer-loop. The three experiments shown are from a strong constraint calculation in primal space (I4D-Var, red curve), a strong constraint calculation in dual space (R4D-Var, blue curve), and a weak constraint calculation in dual space (R4D-Var, black curve). The dashed curve shows the theoretical minimum value of the cost function, Jmin, that will be reached only if all of the prior hypotheses embodied in D are correct

### Number of inner-loops

Fig. 14.3 Log10(J(5z)) (from (14.16)) versus the number of inner-loops for three 4D-Var assimilation calculations for the period 1-4 July, 2000 for the case of 1 outer-loop. The three experiments shown are from a strong constraint calculation in primal space (I4D-Var, red curve), a strong constraint calculation in dual space (R4D-Var, blue curve), and a weak constraint calculation in dual space (R4D-Var, black curve). The dashed curve shows the theoretical minimum value of the cost function, Jmin, that will be reached only if all of the prior hypotheses embodied in D are correct

Fig. 14.3 indicates that as expected from the equivalence of (14.18) and (14.19), both calculations yield the same solution. A comparison of the most likely CCS circulation estimates from 4D-Var in both cases (not shown) confirms that the primal and dual form ultimately yield the same circulation after 75 inner-loops. In both cases, however, the minimum value reached by J(5z) is larger than Jmin indicating that either the prior hypotheses D and R are incorrect, or that the global minimum value of JNL(z) has not been reached.

The influence of different combinations of the number of outer-loops and the number of inner-loops on the final value of JNL(z) in (14.15) that is reached at the end of the 1-4 July, 2000 strong constraint l4D-Var cycle is shown in Table 14.1. Table 14.1 reveals that updating the circulation estimate about which the tangent linear and adjoint models are linearized is clearly beneficial for reducing the value of JNL(z) for the same computational effort (i.e. when the total combined number of outer-loops and inner-loops is the same). In addition, increasing the number of outer-loops yields a value of JNL(z) that is closer to Jmin.

Table 14.1 Values of JNL(z) from (14.15) that are reached at the end of a 4 day strong constraint I4D-Var data assimilation experiment spanning the period 1-4 July, 2000 for different combinations of the number of outer-loops, n, and number of inner-loops, m. In each case, the total number of iterations, n x m, is the same. The value of J . in each case is 0.52 x 104

1 x 100

1.56 x 104

 0.12 0.1 L_ o 0.08 en 0.06 CL 0.04 0.02 0 2 1.5 n LU 1 C/> > CL 0.5 0 01/02/99 01/02/00 01/01/01 01/01/02 01/01/03 01/01/04 12/31/04 01/02/99 01/02/00 01/01/01 01/01/02 01/01/03 01/01/04

12/31/04

Fig. 14.4 Time series of the root mean square (rms) difference between the model and observations for SST and SSH for period Jan 1999-Dec 2004. The blue curve shows the rms differences averaged over 14 day intervals from the model forced only by the prior forcing and subject to the prior boundary conditions. The red curve shows results from a case where primal I4D-Var was applied sequentially over 14 day assimilation windows, and the rms differences are the average over each 14 day window. The control vector was comprised of only the initial conditions in this case. The green curve shows the rms difference associated with the prior circulation estimate xb(t), and is also a 14 day average. Since the prior circulation is the result of a model run started from the most likely circulation estimate at the end of each I4D-Var window, it can also be thought of as a "forecast" of the circulation for the subsequent 14 day period

01/02/99 01/02/00 01/01/01 01/01/02 01/01/03 01/01/04

12/31/04

Fig. 14.4 Time series of the root mean square (rms) difference between the model and observations for SST and SSH for period Jan 1999-Dec 2004. The blue curve shows the rms differences averaged over 14 day intervals from the model forced only by the prior forcing and subject to the prior boundary conditions. The red curve shows results from a case where primal I4D-Var was applied sequentially over 14 day assimilation windows, and the rms differences are the average over each 14 day window. The control vector was comprised of only the initial conditions in this case. The green curve shows the rms difference associated with the prior circulation estimate xb(t), and is also a 14 day average. Since the prior circulation is the result of a model run started from the most likely circulation estimate at the end of each I4D-Var window, it can also be thought of as a "forecast" of the circulation for the subsequent 14 day period

An example of the overall performance of 4D-Var when applied sequentially is shown in Fig. 14.4 which shows the root mean square (rms) difference between ROMS-CCS and the observations for SST and SSH. The result from three calculations are shown corresponding to: (1) no data assimilation, when ROMS-CCS is forced only by the prior forcing fb( t) and subject to the prior boundary conditions bb(t); (2) primal I4D-Var applied sequentially using 14 day assimilation windows; and (3) the prior circulation estimates xb(t) which arise from runs of the model initialized from the most likely circulation estimate at the end of the previous assimilation cycle. The period considered in each case is 1 Jan 1999-31 Dec 2004. Figure 14.4 shows that I4D-Var has a positive impact on the agreement of the model with the observations during each data assimilation cycle, and on agreement with the prior.

### 14.4.4.2 The Control Vector Influence

In the strong constraint examples of 4D-Var presented in Fig. 14.3, the control vector z was comprised of the model initial conditions, x(t0), the surface forcing, f(t), and the boundary conditions, b( t), for the interval t = [t0,tN]. However, it is instructive to explore the relative impact of each component of z on the cost function. To this end, Fig. 14.5 shows the results of a sequence of 4D-Var calculations in which the control vector was successively augmented with different control variables.

x 104

 ■ No assim ■ i.c. only ■ i.c.,wind ■ i.c.,wind,Q □ i.c.,wind,Q,(E-P) □ i.c.,wind,Q,(E-P),b.c. ■ i.c.,wind,Q,(E-P),b.c., weak

Experiment

Fig. 14.5 The cost function J(dz) in (14.16) after 1 outer-loop and 75 inner-loops from several dual R4D-Var experiments for the 4 day period 1-4 July, 2000. Expt. A: no data assimilation; Expt. B: dz comprised of the initial condition increments dx(t0) only; Expt. C: dz comprised of only dx(t0) and the surface wind stress components of df(t); Expt. D: dz comprised of only dx(t0) and the surface wind stress and heat flux components of df(t); Expt. E: dz comprised of dx(t0) and all components of df(t); Expt. F: dz comprised of dx(t0), df(t), and db(t), Expt. G: same as F except using the weak constraint as described in Sect. 14.4.4.3. The dashed line indicates the theoretical minimum value of the cost function, J

Figure 14.5 reveals that as control vector increment Sz is augmented with more components of the surface forcing and boundary conditions, the minimum value reached by the cost function J(Sz) progressively decreases. This is because the length of Sz is increasing at the same time meaning that there are more degrees of freedom available to 4D-Var with which to fit the observations. It is important to realize, however, that the 4D-Var process is not linear, so if the order of the calculations in Fig. 14.5 is changed, the same progression of changes in J(Sz) will not necessarily be obtained. Nonetheless, Fig. 14.5 does suggest that accounting for the errors and uncertainties in the initial conditions has by far the largest impact on the efficacy of the circulation estimate measured in terms of J(Sz). Significant further reductions in J(Sz) are afforded when uncertainties in surface forcing are accounted for also. However, Fig. 14.5 suggests that uncertainties in the boundary conditions have the least impact on J(Sz).

### 14.4.4.3 Weak Constraint 4D-Var

Model errors are a significant source of uncertainty in model-derived circulation estimates, and it is important to account for these uncertainties during 4D-Var. However, quantifying and identifying the sources of model error is one of the greatest challenges in ocean data assimilation. In ROMS-CCS, we have made some progress towards identifying some of the most significant impacts of model error, but our understanding is far from complete.

A dominant feature of the CCS circulation is the occurrence of coastal upwelling in the spring and summer along much of the California, Oregon and Washington coasts which is driven by equatorward alongshore winds. ROMS-CCS simulates the seasonal cycle of upwelling very well (Veneziani et al. 2009) although in the absence of data assimilation, comparisons of the model with observations indicates that the model SST is biased compared to the observed temperatures during the peak of the upwelling season. Independent investigations of the quality of the surface wind forcing by Doyle et al. (2009) have revealed that the COAMPS priors, fb(0, used to drive ROMS-CCS agree well with wind observations from satellite scatterometers and ocean buoys. Therefore, it seems likely that the model bias in ROMS-CCS SST is associated with model errors rather than errors in the surface forcing. Further evidence for this hypothesis is provided by the work of Broquet et al. (2009b, 2011) using strong constraint 4D-Var. They found that data assimilation leads to a reduction in strength of the upwelling favorable alongshore surface winds during spring and summer, and a general degradation of the wind prior when compared to satellite-derived surface wind observations. An important conclusion of their work, therefore, is that in the absence of corrections for model error, data assimilation may yield undesirable and non-physical corrections to potentially all elements of the control vector in an attempt to minimize the cost function.

Based on the findings of Broquet et al. (2009b, 2011) attempts are currently underway to account for model error in ROMS-CCS during 4D-Var using the weak constraint. The impact of imposing a weak constraint during R4D-Var on the cost function J(Sz) is also shown in Fig. 14.3. Based on the known bias of the ROMS-CCS SST in spring and summer, it was assumed that the model error affecting SST

is present mainly in the model temperature equation, and confined to within a short distance of the coast. To account for such errors during 4D-Var, then using (14.14) we let x(ti) = M(ti,ti-i)(xt(ti-i), f (ti),b(ti)) + eq at each model time step where xt, ft, and bt are the true state, forcing, and boundary conditions, and eq represents the model error. A model error covariance matrix of the form Q = [eq, e^) = WB was assumed, where B is the initial condition prior error covariance matrix, and W is a diagonal rescaling matrix. All elements of W are zero except those corresponding to ocean temperature grid points within 300 km of the North American coast. The non-zero elements of W are of the form a(1-d/300), where a = 0.05 is a variance scaling factor, and d is the distance measured from the coast in km. The choice for a represents a standard deviation for the prior model error £q in temperature that is ~22% of that of the error in the initial condition prior.

Figure 14.3 shows that the weak constraint R4D-Var cost function J(bz) asymptotes to a lower value compared to the strong constraint case, and is closer to the theoretical minimum value Jmin. Additional experiments have revealed that the minimum value of J(bz) reached during weak constraint R4D-Var decreases with increasing a. However, for a > 0.05 the resulting circulation estimates possess features that are very unphysical. Nonetheless, these preliminary attempts at accounting for uncertainties due to model during 4D-Var are encouraging.

Figure 14.5 also shows the impact of a weak constraint on J(bz) in relation to the sequence of experiments described in Sect. 14.4.4.2 where the control vector is successively augmented. Experiment G is the same as that shown in Fig. 14.3, but only the minimum value reached by J(bz) is shown. Figure 14.5 indicates that augmenting the increment control vector bz with n( t) leads to a further reduction in J(bz), and J(bz) is closer still to Jn

### 14.4.5 4D-Var Diagnostics

Having successfully computed the most likely circulation estimate using 4D-Var, there are a number of diagnostic calculations that can be performed using the 4D-Var output that are of considerable interest. Two examples are presented here that provide information about the accuracy of the resulting circulation in the form of the posterior error covariance, and information about the impact of each individual observation on different physical aspects of the circulation.

The calculations presented below are greatly facilitated by the Lanczos formulation of the conjugate gradient algorithm that is used to minimize J(bz) in the ROMS 4D-Var algorithms (Fisher and Courtier 1995). Specifically, the primal form of the gain matrix (14.18) can be expressed as K = VpTp ^T GTR-1, where each column of the matrix Vp is a primal Lanczos vector (Golub and van Loan 1989), and Tp is a tridiagonal matrix. Each inner-loop yields one additional member of the Lanczos sequence which form an orthonormal basis, and after m inner-loops the matrix Vp has dimension Nm x m. Similarly, the dual form of the gain matrix (14.19) can be expressed as K = DGTVdTp1 V where Vd is a matrix of dual Lanczos vectors of dimension N , x m.

14.4.5.1 Posterior Error

If zt represents the control vector describing the true state of the ocean, then the error covariance matrix of the posterior estimate za=zb + ¿za is given by Ea=E((za- zt) (za- zt)T) where E is the expectation operator. It is easy to show that in the case where the observation errors and prior errors are uncorrelated, the posterior error covariance matrix can be expressed as:

where K is the gain matrix of Sect. 14.3.2.

Exercise 6: Using the definitions Ea=E((za-zt)(za- zt)T) and D=E((zb- zt)(zb- zt)T) and the dual form of K in (14.19), derive Eq. (14.24) for the posterior error covariance matrix assuming uncorrelated prior errors and observation errors.

Figure 14.6 shows an example of the reduction in the uncertainty in SST and subsurface temperature as a result of data assimilation using R4D-Var. In this case, ROMS-CCS with 10km horizontal resolution and 42 vertical levels was used. The diagonal elements of D and Ea represent the prior and posterior error variance respectively of each control variable. Therefore, the diagonal of the difference matrix A=(Ea- D) represents the change in error variance of the prior due to 4D-Var. Negative values of the diagonal of A represent posterior grid point variables that are more certain than the prior, while zero values indicate that the posterior and prior are equally uncertain. Figure 14.6 shows the diagonal elements of A corresponding to SST and 75 m temperature, and indicates that the posterior upper ocean tem- ,0.00 -0.03 -0.06 -0.09 -0.12 -0.15 -0.18 -0.21 -0.24 -0.27 Fig. 14.6 Posterior variance-prior variance for SST (left) and 75 m temperature (right) for a typical R4D-Var calculation in ROMS-CCS with 10 km resolution

perature is more certain than that of the prior over large areas of the ROMS-CCS domain, particularly in the region of high eddy kinetic energy along the California central coast identified in satellite observations by Kelly et al. (1998). However, a degree of caution must be exercised when interpreting Ea using the reduced rank approximation of K based on the Lanczos vectors described above. Ea in (14.24) describes the expected posterior error covariance only in the case where the priors D and R are correct, and when 4D-Var is run to complete convergence. In general, neither of these conditions will be satisfied, and since the number of inner-loops m is typically « Nobs, the Lanczos vectors Vd may span only a small portion of observation space. This results in an over estimate of the posterior error variances as demonstrated by Moore et al. (2011b).

### 14.4.5.2 Observation Impacts

It is of considerable interest to know which observations or which types of observations exert the greatest influence on particular physical aspects of the most likely circulation estimate derived from 4D-Var. For example, Fig. 14.7 shows a vertical section of the time average of the alongshore velocity over the upper 500 m of the water column along 37N for the period 1 Jan 1999-31 Dec 2004 from two ROMS-CCS calculations: one in which no data were assimilated, and one in which data were assimilated sequentially every 14 days using strong constraint I4D-Var. In the latter case, the increment control vector is comprised only of dx(t0) as described in Broquet et al. (2009a). Clearly the assimilation of observations has a significant impact on the structure of the CCS at this latitude.

In order to quantify the impact of each individual observation on the circulation, consider the transport at time t across the section shown in Fig. 14.7 for the posterior circulation of a typical 4D-Var cycle, namely Ia(t) = hTxa(t), where h is a vector with zero elements everywhere except for those elements that correspond to velocity grid points that contribute to the transport normal to the section. The non-zero elements take the form of the transformation coefficients required to rotate the total velocity in the alongshore direction, and the appropriate area elements for each grid cell in the vertical. However, recall that the circulation estimate xa(t)=xb(t) + dxa(t), in which case Ia(t) = hTxb(t) + hTdxa(t) =Ib(t) + hTdxa(t). Therefore, the difference in transport AI(t) = Ia(t)-Ib(t) between the posterior and the prior circulation estimates is given by AI(t) = hTdxa(t). However, recall that to first-order dxa(t) = M(t,t0)Kd where M(t,t0) is the tangent linear model. The difference in transport between the posterior and the prior can then be expressed as:

where MT(t0,t) is the adjoint model. For each data assimilation cycle, (14.25) can therefore be used to compute the contribution of each observation, represented by the individual elements of d, to the change AI(t) in the prior transport associated with 4D-Var. a -0.1 -0.05 0 0.05 0.1 Fig. 14.7 Summer time average (July-September) vertical sections of the alongshore component of velocity from ROMS-CCS during the period 1 Jan 2000-31 Dec 2004 for the case a where there is no data assimilation, and b where data are assimilated sequentially every 14 days using strong constraint I4D-Var. (From Broquet et al. 2009a)

Fig. 14.7 Summer time average (July-September) vertical sections of the alongshore component of velocity from ROMS-CCS during the period 1 Jan 2000-31 Dec 2004 for the case a where there is no data assimilation, and b where data are assimilated sequentially every 14 days using strong constraint I4D-Var. (From Broquet et al. 2009a)

Figure 14.8 shows an example calculation from a 7 day strong constraint R4D-Var calculation for the period 5-11 April, 2003. Figure 14.8b indicates that the change in the prior estimate transport along 37N on 11 April as a result of assimilating ~1.5 x 104 observations is ~0.75Sv. A comparison of Fig. 14.8a, b, however, reveals that even though satellite 