## Abstract

A computationally efficient relocatable system for generalized inverse (GI) modeling of barotropic ocean tides is described. The GI penalty functional is minimized using a representer method, which requires repeated solution of the forward and adjoint linearized shallow water equations (SWEs). To make representer computations efficient, the SWEs are solved in the frequency domain by factoring the coefficient matrix for a finite-difference discretization of the second-order wave equation in elevation. Once this matrix is factored representers can be calculated rapidly. By retaining the first-order SWE system (defined in terms of both elevations and currents) in the definition of the discretized GI penalty functional, complete generality in the choice of dynamical error covariances is retained. This allows rational assumptions about errors in the SWE, with soft momentum balance constraints (e.g., to account for inaccurate parameterization of dissipation), but holds mass conservation constraints. While the dynamical calculations involve elevations alone, depth-averaged currents can be directly assimilated into the tidal model with this approach. The efficient representer calculation forms the basis for the Oregon State University (OSU) Tidal Inversion Software (OTIS). OTIS includes software for generating grids, prior model covariances, and boundary conditions; for time stepping the nonlinear shallow water equations to generate a first guess or prior solution; for preliminary processing of TOPEX/Poseidon altimeter data; for solution of the GI problem; and for computation of posterior error bars. Approximate GI solution methods, based on using a reduced set of representers, allow very large datasets to be inverted. OTIS regional and local GI tidal modeling (with grids containing up to 10^{5} nodes) require only a few hours on a common desktop workstation. Use of OTIS is illustrated by developing a new regional-scale (1/6°) model of tides in the Indonesian Seas.

## 1. Introduction

A large fraction of the variance in many oceanographic variables is due to tides. As a result, for many practical applications in the marine environment (e.g., shipping and modeling of pollution dispersal) accurate maps of tidal currents or elevations are often indispensable. Knowledge of tides is of value in many scientific studies as well. Here the tides are often seen as a nuisance that must be removed from the data prior to studies of longer period oceanographic variations. For a fixed observation site, such as a semipermanent mooring, tides can be removed by filtering of time series. With data taken from moving platforms (e.g., ships, Lagrangian drifters, satellites) spatial and temporal variations are aliased, rendering a simple filtering scheme untenable, and again requiring accurate maps of tidal fields. Tides are of course not just a source of noise. For example, studies of ocean microstructure have shown a correlation of turbulent dissipation rates with tidal cycles (e.g., Polzin 1997; Ledwell et al. 2000), and there is increasing evidence that tides may provide a significant source of energy for mixing the deep ocean (Munk and Wunsch 1998; Egbert and Ray 2000). For studies of this sort a detailed knowledge of the tidal fields is again desirable.

The basic equations of tidal dynamics (e.g., Hendershott 1977) are comparatively simple and have been understood since the time of Laplace. However, there are a number of complications that still make accurate modeling of even the barotropic tides a challenging problem in practice. Difficulties include the need for accurate open boundary conditions and bottom topography, the need for approximate parameterizations of dissipation in the tidal equations, solid earth effects, and the effects of ocean stratification on the barotropic tides, which may be difficult to account for without full 3D modeling of baroclinic tidal currents. Examples of the last complication include extraction of barotropic energy due to scattering over rough bottom topography into baroclinic waves (Bell 1975; Baines 1982) and effects of stratification in deep water on depth-independent tidal currents on the continental shelf (Cummins et al. 2000).

Data assimilation methods offer a way around some of these difficulties. By requiring the tidal fields to fit estimates of tidal constants at even a limited number of locations in the model domain, accuracy of tidal predictions can be significantly improved. Early applications of some sort of data assimilation scheme include the “hydrodynamic interpolation” scheme of Schwiderski (1980) and early discussions of data assimilation by Bennett and McIntosh (1982) and McIntosh and Bennett (1984). Recent reviews of tidal data assimilation are given by Egbert and Bennett (1996) and Kivman (1997). Our focus here is on the generalized inversion (GI) scheme of Egbert et al. (1994), hereinafter EBF, which was developed for assimilating TOPEX/Poseidon (T/P) altimeter data into a global barotropic tidal model. The approach used involves minimizing a weighted sum of squared misfits to the linearized shallow water equations and the data. To accomplish this, EBF used a “representer” approach (Bennett 1992). This scheme, which we briefly review below, has some advantages over simple schemes (see Egbert and Bennet 1996; Egbert 1997; Kivman 1997) but is very computationally intensive, since the dynamical equations must be solved essentially twice (once forward in time and once backward) for each observation.

For very large datasets such as the T/P altimetry, some modifications of the computational approach are thus required. EBF solved the global tidal inverse problem using a greatly reduced set of representers. Egbert and Bennett (1996) and Egbert (1997) developed an iterative solution method based on minimization of the penalty functional in the data space using conjugate gradients. This approach still uses a representer formulation, but for large datasets it requires many fewer solutions of the dynamical equations. Even with these refinements to the representer method, inverse modeling of the large T/P dataset on a ⅔° nearly global grid required a major computational effort. Several months of CPU time on a small supercomputer (a 32-node CM-5) was required for the solution (Egbert 1997).

Here we describe a much more efficient numerical implementation of a GI scheme for barotropic tidal modeling. The basic idea is quite simple: greatly improve the efficiency of the numerical modeling scheme. EBF solved the linear shallow water equations by explicit time stepping of the momentum and continuity equations with a periodic forcing, followed by harmonic analysis of the steady-state solution. Numerical stability requires the time step in such a scheme to be very short compared to the tidal period, so long run times are required. As resolution is increased (e.g., for regional scale instead of global modeling) even shorter time steps and longer run times are required. Here we solve the equations in the frequency domain by direct factorization of the coefficient matrix for a finite-difference approximation of the dynamical equations. The advantage of this approach for the inverse problem is clear: once the coefficient matrix has been factored, repeated solution of both forward and backward equations (as required by the representer approach) is extremely fast.

The major difficulty with this approach is the very large size of the factored coefficient matrix for numerical grids of reasonable size. Eliminating velocities from the shallow water equations to yield a second-order elliptic equation for the tidal elevation alone reduces the number of unknowns by a factor of 3, and thus matrix size by a factor of 9. A similar approach has been used previously in the finite-element inverse scheme described by Lyard (1999). A key feature of the approach developed here is that we retain the first-order system of shallow water equations (defined in terms of elevations and currents) in the definition of the GI penalty function and only eliminate variables for the representer calculation. With this formulation, complete generality in definition of the dynamical error covariance is retained, for example, allowing errors in the momentum equations, but requiring exact conservation of mass (see EBF). Thus currents can be computed from elevations without assuming that the depth-integrated momentum balance equations (which depend on bathymetry and parameterization for dissipation and solid earth effects) are exact; and current data can be assimilated, even though all dynamical calculations are in terms of elevations alone. This general approach should be applicable to other oceanographic data assimilation problems where it may be convenient to eliminate some variables for dynamical calculations.

We have combined this efficient representer calculation scheme with software for generating grids, prior model covariances, boundary conditions, and reduced altimetry datasets into a relocatable package of programs for tidal data inversion. With this package, which we refer to as OTIS (for OSU Tidal Inversion Software), regional-scale tidal inverse problems (with grids containing up to 10^{5} nodes) can be solved in a few hours on a common desktop workstation with a reasonable amount (say 500 MB) of RAM.

Our primary goals in this paper are to present the efficient representer calculation method and to provide an overview of the relocatable tidal inversion software package OTIS. In section 2 we review the tidal inverse problem and the representer approach. In section 3 we describe our approach for efficient calculation of representers, both for elevation and current data. In section 4 we provide some further details on OTIS, including methods for dealing with the very large T/P altimeter dataset and calculation of posterior errors. Finally, in section 5 we provide an illustrative application of OTIS, inverting for barotropic tides in the seas surrounding Indonesia.

The philosophy and theoretical underpinnings of our inversion approach are not given extensive discussion here. The reader is referred to EBF, Egbert and Bennett (1996) and Egbert (1997) for this perspective.

## 2. Generalized inversion of tidal data

Our goal is to find tidal fields **u**, which are consistent both with the hydrodynamic equations, which we denote formally as

and with a *K*-dimensional vector **d** of tidal data

In (2) 𝗟 = [*L*_{1}, … , *L*_{K}] corresponds to the *K* measurement functionals, which relate the observed data to the unknown tidal state **u**. Due to measurement errors and inadequacies in the necessarily approximate dynamical equations, there will be in general no **u** that exactly satisfies both equations. With the GI approach we compromise between (1) and (2) by minimizing the quadratic penalty functional

Here **Σ**_{e} and **Σ**_{f} are the covariances for the data and dynamical errors, which express our a priori beliefs about the magnitude and correlation structure of errors in the data, and in the assumed dynamical equations.

Specification of the dynamical error covariance **Σ**_{f} is a critical part of the inverse problem, but is not the focus of our discussion here. For all examples discussed (and for OTIS) we use a covariance of the general form discussed in EBF. This allows for spatially varying dynamical error amplitudes, with a constant decorrelation length scale. Amplitudes are estimated as described in EBF based on an analysis of the dynamical equations, using a preliminary global tidal model to estimate typical amplitudes of the tidal fields. Further details will be given for the specific examples of section 5. As in EBF, we assume a simple diagonal form for the data error covariance **Σ**_{e}.

### a. The representer approach

If the dynamical equations of (1) are linear, the representer approach (Bennett 1992) can be used to minimize (3). The minimizer of (3) can be written as

where **u**_{0} = 𝗦^{−1}**f**_{0} is the exact solution of (1) and the functions **r**_{k} are the representers (Yosida 1980) of the data functionals defined by *L*_{k}, *k* = 1, *K.* Representers can be calculated by first solving the adjoint of the dynamical equation

where Δ_{k} is the averaging kernel for the data functional *L*_{k} (an impulse at **x**_{k} for a point observation at this location), and then solving the forward equation

Note that the forcing for (6) is the solution to (5) smoothed by convolution with the dynamical error covariance **Σ**_{f}. The representer coefficients ** β** = (

*β*

_{k},

*k*= 1,

*K*), are found by solving the

*K*×

*K*system of equations

where 𝗥 is the representer matrix with elements

In the most straightforward application of the representer approach, (5) and (6) are solved *K* times each for the functions **r**_{k}, 𝗥 is constructed using (8), (7) is solved for ** β**, and finally the solution

**û**is calculated via (4). As discussed in EBF, Egbert and Bennett (1996), and Egbert (1997), there are several variants on this approach that approximate the full solution and are considerably more efficient for large datasets. We consider some practical aspects of these variants further in section 4.

### b. Dynamical equations

We assume essentially linear shallow water dynamics of the form

Here *ζ* is the elevation of the sea surface; **U** is the volume transport vector, equal to velocity times water depth *H*; *f* is the Coriolis parameter; **ẑ** is oriented to the local vertical; and **F** is the frictional or dissipative stress. The astronomical tide-generating force with allowance for the earth's body tide (Hendershott 1977) is denoted by **f**_{0}.

Tidal loading and self-attraction (Hendershott 1972; Ray 1998) are accounted for by the term *ζ*_{SAL}. EBF used the crude approximation *ζ*_{SAL} = 0.1*ζ* (e.g., Schwiderski 1980). In fact *ζ*_{SAL} should be computed by convolution of *ζ* with the Green's function for loading and self-attraction (Farrell 1972; Ray 1998). Since this convolution smoothes out small-scale features, and since large-scale tidal elevations are now well determined over most of the earth, *ζ*_{SAL} is in fact now reasonably well known, even where local details of tidal elevations and currents remain uncertain. We thus move *gH*∇*ζ*_{SAL} to the right-hand side of (9), so that the total specified forcing becomes **f** = **f**_{0} + *gH*∇*ζ*_{SAL} with *ζ*_{SAL} computed from a global tidal model. Note that since we ultimately do not fit the dynamics of (9) exactly, small errors due to this approximate treatment of tidal loading and self-attraction can be corrected by the data. Boundary conditions for (9) are no flow across the coast, and specifications of either elevations or the normal component of volume transports on any open boundaries. As with the dynamical equations, we allow for errors in all of the boundary conditions.

EBF assumed a simple linear parameterization of the dissipation term **F** = (*r*/*H*)**U**, where *r* is the bottom friction coefficient. A more physically reasonable parameterization of dissipation due to bottom boundary layer drag is quadratic in velocity (e.g., Pingree 1983)

where **v** is the total velocity vector (in particular, including all tidal constituents), and the value of the nondimensional parameter *c*_{D} is approximately 0.003. Parameterizing the dissipation as in (11) makes the shallow water equations nonlinear, so that the representer theory for minimizing (3) is not strictly applicable. We use a simple linearization approach. The prior solution **u**_{0} is computed for all major tidal constituents simultaneously by time stepping the shallow water equations with the quadratic dissipation of (11). We also include the inertial terms (**U** · ∇**v**) in (9), the nonlinearity in (10) (i.e., *ζ***U**), and horizontal eddy viscosity in this computation. For most applications to regional-scale models (with grid resolutions down to 10 km or so) including these additional terms appears to make little difference to the final tidal solutions, and will thus not be discussed further here.

With the prior solution computed we linearize (11). The simplest approach is to calculate the time-averaged tidal velocity (including all constituents) to yield a spatially varying linear drag coefficient *κ* = (*c*_{D}‖**v**‖/*H*). Several refinements are possible. The time domain equations can be linearized around the solution **u**_{0}, Fourier transformed, and simplified by retaining only the significant tidal frequencies. In the resulting linearized operator the dissipative stress takes the form **F** = 𝗞**U** where 𝗞 is a spatially varying drag tensor. There is also a weak coupling in the linearized equations between constituents. Other approaches to linearizing the dissipation using more rigorous asymptotic analysis are discussed by Le Provost and Poncet (1978). We will not pursue these refinements here, since all terms are derived from an approximate (and to some degree ad hoc) parameterization of **F**. Again we do not fit the dynamics exactly, and errors due to our approximate treatment of dissipation can be compensated by the data. With the simple linearization used here, the dissipative terms in the dynamical equations are at least of the proper order of magnitude (in contrast to the simpler constant linear drag coefficient used by EBF).

With dissipation linearized, the dynamical equations are linear, and can be transformed to the frequency domain with individual tidal constituents (e.g., *M*_{2}, *S*_{2}, *K*_{1}, *O*_{1}) decoupled. For our initial discussion of the representer calculation we thus take the tidal state vector **u** = (*U, **V, **ζ*) to be volume transports and elevations for a single tidal constituent of frequency *ω.* We take 𝗦 to be the linearized frequency domain shallow water equations, and we assume that the data are provided as harmonic constants of either elevations or currents at a series of locations **x**_{k}. In this case the forcing for the representer calculation (8) (i.e., *L*_{k}) corresponds to an impulsive forcing at **x**_{k}. As discussed in EBF, for direct inversion of time domain data (as might be appropriate for altimeter data with short time series, or shipborne acoustic Doppler profiler measurements of currents where harmonic analysis at fixed locations may be impossible) the observation functionals couple all constituents together. Representers for this more complicated case can be constructed from the single constituent representers that we focus on here (EBF). We summarize the implementation of this complication into OTIS in section 4.

## 3. Efficient calculation of representers

The bulk of the calculation required for the GI approach of EBF lies in computing the representers. Here we describe an efficient scheme for the repeated solution of the shallow water equations required for this calculation. With the notation and approximations of the previous section, the frequency domain equations may be written as

where

For now it is useful to allow for arbitrary forcing and inhomogeneous boundary conditions of both the momentum and continuity equations. Appropriate forcing and boundary conditions will be discussed more explicitly below. Assuming **Ω** is invertible at all locations (this will be true providing *κ* ≠ 0 or *ω* ≠ *f*), we can rewrite (12) as

and combine this with (13) to get a second-order equation in *ζ*:

### a. Numerical implementation

To describe our numerical implementation of the representer calculation, we introduce the discrete operators 𝗚,𝗗,𝗖, and𝗔, which correspond, respectively, to

For the discrete approximation we use a C-grid (Fig. 1). On this staggered grid 𝗚 maps naturally from a pair of adjacent *ζ* nodes onto the intermediate *U* or *V* node. Similarly, 𝗗 maps onto *ζ* nodes from the four surrounding *U* and *V* nodes. The operator 𝗖, which maps from *U* and *V* nodes back to *U* and *V* nodes is not so naturally defined on a C grid. To implement this operator on a *U*(*V*) node requires first averaging over the four adjacent *V*(*U*) nodes. Details are given in appendix A. The operator 𝗔 defined in terms of 𝗚, 𝗗, and 𝗖 (and multiplication by *iω*), maps from *ζ* nodes to *ζ* nodes and represents the second-order operator of the frequency domain wave equation.

In terms of these discrete operators the full solution **u** = (*ζ, **U*) to (1) can be expressed formally as

Arbitrary inhomogeneous boundary conditions at the coast and open boundaries are readily implemented through appropriate definitio of **𝗚,𝗗,𝗖,𝗔**, **f**_{ζ}, and **f**_{U} on the boundary nodes. Details are given in appendix A.

The partitioned matrix on the right side of (18) is a discretization of the operator 𝗦^{−1} needed to solve (6). From (18) it is simple to derive the form for the discrete inverse adjoint operator required to find the solution to (5). For a general forcing Δ = (Δ_{ζ}, Δ_{U}), this takes the form:

where the superscript asterisk denotes the conjugate transpose of the matrix. To solve (5) and (6), we thus proceed as follows. First the coefficient matrix 𝗔 for the wave equation is formed and factored into lower- and upper-triangular matrices as 𝗔 = 𝗟𝗨. The full matrix has a bandwidth of (*m* + 1) and, allowing for full pivoting, requires (3*m* + 4)*n**m* complex words to store (Anderson et al. 1992). For a global model domain (periodic in longitude) the matrix is no longer strictly banded, but the solution algorithm is easily modified to allow for this case (see appendix B).

Since 𝗔* = 𝗨*𝗟*, multiplication byboth 𝗔^{−1} and 𝗔^{−1*} [required to use (18) and (19), respectively, to solve (5) and (6)] can be accomplished rapidly by solving the appropriate triangular systems. All of the other discrete operators appearing in (18) and (19) are very sparse (and local) and hence require minimal computation. With 𝗔 factored and stored in RAM, we can rapidly compute a large number of representers, corresponding to observations at different locations.

The above discussion is completely general, allowing for solution of the shallow water equations and adjoint with arbitrary forcing. It is instructive to explicitly consider the simplifications that result for some more specialized (but highly relevant) cases.

### b. Representers for elevation data

First we consider harmonic constants for tidal elevation at locations **x**_{k}, *k* = 1, *K.* In this case the averaging kernel for data functional *L*_{k} takes the form of an impulse in the *ζ* component at **x**_{k}. In (19) we thus have Δ_{U} = 0. Further simplification in the representer calculation results from assuming, as in EBF, that the continuity equation (13) is exact. With the dynamical error covariance **Σ**_{f} partitioned between *ζ* and *U* − *V* nodes as in (18), this implies

that is, **Σ**_{f} acts only on *U* and *V* nodes. We can thus formally write the elevation component of the representer as

where Δ_{k} is a unit magnitude impulsive forcing at **x**_{k}. The transport component of the representer **r**_{U} is also readily calculated, but if only elevation data is inverted, computation of the representer matrix 𝗥 needed to solve (7) for the coefficients ** β** requires only

**r**

_{ζ}.

Explicitly the steps implied by (20) are the following:

Using the LU factorization of 𝗔, solve the conjugate transpose system with impulsive forcing at

**x**_{k}: 𝗔**α*_{ζ}= Δ_{k}. The solution is defined on the*ζ*nodes.Apply the local operators 𝗗* and then 𝗖*. These are the transpose of 𝗗 and 𝗖, respectively, and map

*α*_{ζ}to*α*_{U}, defined on the*U*/*V*nodes.Apply the covariance smoother

**Σ**_{f}to*α*_{U}.Aply the local operators 𝗖 and then 𝗗. This maps the smoothed forcing of the momentum equations (

**Σ**_{f}*α*_{U}) to a forcing of the wave equation defined on the*ζ*nodes.Using this forcing, solve the wave equation for

**r**_{ζ}.

One issue we have not yet discussed explicitly is the specification of boundary conditions for the representer calculation. As shown in EBF, boundary conditions for (5) are homogeneous (unless the data involve a measurement on the boundary), but inhomogeneous boundary conditions are required for (6). EBF give expressions for the appropriate boundary conditions in terms of the adjoint solution and boundary condition error covariances. If the discrete operators 𝗚,𝗗,𝗖, and𝗔 are set up as described in appendix A to allow for arbitrary elevations on boundary *ζ* nodes, and arbitrary normal volume transport on boundary *U* and *V* nodes, the appropriate boundary conditions will automatically be computed at steps (i) (for elevation boundary nodes) and (ii) (for *U* and *V* boundary nodes). We discuss this further in appendix A.

As in EBF, separate covariances are used for smoothing and scaling coastal boundary conditions, open boundary conditions, and the forcing for (6). Errors in the coastal boundary conditions are assumed to arise solely from approximation of the coast on the numerical grid, and are thus not smoothed. Convolution with appropriate boundary condition covariances is implicit in step (iii) above, which is identical to the scheme described in EBF. Since the model covariance **Σ**_{f} is symmetric, it is clear from (20) that the representer matrix must be Hermitian (to the precision of numerical round-off error). This is an advantage of the direct calculation scheme used here. Numerical implementation of the forward and adjoint problems, including all boundary conditions, is guaranteed to be consistent.

### c. Calculation of the inverse solution: Currents

With the elevation components of the representers calculated, we may construct 𝗥 as in (8) (i.e., *R*_{jk} is computed by evaluating **r**_{k} at **x**_{j}) and solve (7) for ** β**. Rather than explicitly forming the sum of representers to calculate the GI solution as (4), we redo the equivalent of one representer calculation (8) with the single impulsive forcing Δ

_{k}replaced by the combination Σ

_{k }

*β*

_{k}Δ

_{k}(see EBF).

Repeating steps (i)–(iii), we obtain

This is exactly the forcing error required to bring the inverse solution **û** = **u**_{0} + *δ***u** into agreement with the data. Note that with the model covariance of EBF, *δ***f** corresponds to forcing of the momentum equations, together with coastal and open boundary conditions. The elevation component of *δ***u** is calculated from *δ***f** as in steps (iv) and (v) of the representer calculation as

From (18) the corresponding volume transports are readily seen to be

It is important to note that the last step requires knowing the forcing error for the momentum equations. The correction to the volume transport component *δ***U** cannot be recovered from *δζ* alone, because the momentum equations, which relate *δ***U** to the gradient of *δζ,* do not hold exactly.

Near the critical latitude where *ω* = *f* the 2 × 2 matrices **Ω**^{−1} will be nearly singular if the drag coefficient *κ* is small. This can result in noisy estimates of *δ***U**, especially when calculations are done in single precision. In some cases, numerical round-off error can result in transports *δ***U**, which are not consistent with the continuity equation. To stabilize calculation of *δ***U**, we thus replace the direct calculation of (23) with a weighed least squares fit of *δ***U** to both (23) and ∇·*δ***U** = −*iωδζ.* A scheme similar to that described by Ray (2001) is used.

### d. Representers for velocity data

For a component of the depth-averaged velocity measured at **x**_{k} the forcing for the adjoint equation (19) is of the form Δ = (0, Δ_{U}), where Δ_{U} is a unit magnitude impulse, divided by water depth *H* (to account for our formulation in terms of velocity transports). Again assuming that the continuity equation is exact, (18) and (19) imply that the representer can be calculated as

In this case it is most convenient to calculate the full representer (elevation and transport components), although if only current measurements are to be inverted, only the transports need be saved for calculation of the representer matrix.

The velocity representer calculation is somewhat more complicated, but the basic steps are as outlined above. There are two differences to note in (24). First, the impulse Δ_{U} must be converted to a forcing for the adjoint wave equation 𝗚*𝗖*Δ_{U}. Second, there is an additional term 𝗖*Δ_{U} (calculated directly from the impulse forcing) that must be added to *α*_{U} after solving the adjoint equation, but before smoothing by **Σ**_{f}.

## 4. Solution of large tidal inverse problems with OTIS

A literal application of the representer approach outlined in section 2 requires computation of *K* representers (one for each data location), followed by solution of the *K* × *K* system of equations (7), and formation of the sum of (4). For inversion of altimeter data, where *K* can easily exceed 10^{5}, a direct application of this approach is impractical. In this section we consider approximations to the representer approach, which are practical for large datasets. We also discuss some other ancillary aspects of the tidal inverse calculations, including initial data processing steps for altimeter data, treatment of time domain data, and calculation of posterior errors. For the most part we use variants on methods that have been described previously in some form in EBF or Egbert and Bennett (1996), so we present only a summary of the implementation of these calculations in OTIS.

We initially consider inversion of harmonic constants of tidal elevations and/or currents for a single constituent. We then describe calculation of harmonic constants for altimeter data, and then the extension to directly fitting time domain data where multiple constituents are coupled and must be treated simultaneously.

### a. The reduced basis approach

Representers for measurements of tidal elevations or currents at nearby locations will typically be very similar. More generally, the full set of *K* representers, which form a basis for the minimizer of (3), are highly redundant. It is thus typically possible to form a very good approximation to the minimizer of (3), using only a subset of the representers. With a reduced basis approach (Parker and Shure 1982; EBF), we calculate a relatively small subset of the representers and then minimize the penalty functional over linear combinations of these. Note that although only a subset of representers are computed and used as basis functions for the inverse solution, all available data are still fit (Egbert and Bennett 1996).

Let *K* be the number of locations where harmonic constants are available. We choose a subset *N* of these data locations, and calculate representers for this data subset. Let 𝗥 be the square *N* × *N* matrix defined in (8), with elements corresponding to the *N* calculated representers evaluated at the corresponding *N* data locations. Next, let 𝗣 be the rectangular *K* × *N* matrix with elements *P*_{kj} = *L*_{k}**r**_{j}, that is, the calculated representers evaluated at all *K* of the data locations. We seek the minimizer of (3) among all linear combinations of the *N* calculated representers (and the prior model)

For **u** defined by (25) we clearly have

Also, EBF show

so for the subspace of solutions defined by (25) the penalty functional (3) can be written in terms of the representer coefficients ** β** as

where **d**′ = **d** − 𝗟**u**_{0}. Note that in the limiting case where all representers are calculated 𝗣 = 𝗥, and the coefficients obtained by minimizing (26) yield the global minimum of (3).

To minimize (26) over the reduced coefficient vector ** β**, we compute the singular value decomposition (SVD) of

where 𝗘 = 𝗥^{−1/2}, 𝗪, and 𝗤 are, respectively, *K* × *N* and *N* × *N* matrices with orthonormal columns, and 𝗦 is the diagonal matrix of real positive singular values *s*_{i}. It is simple to show (e.g., Parker 1994) that the minimizer of (26) can then be written as

where *υ* = 1. To calculate 𝗘 = 𝗥^{−1/2} we use an eigenvector decomposition of the representer matrix 𝗥 = 𝗩**Ψ**𝗩*, so that 𝗥^{−1/2} = 𝗩**Ψ**^{−1/2}, where 𝗩 is orthonormal and **Ψ** is the diagonal matrix of positive real eigenvalues. To ensure stability in calculations with 𝗘 we eliminate eigenvectors corresponding to small eigenvalues *ψ*_{i}. This truncation of representers based on the eigenvector decomposition of 𝗥 can also be used to reduce the size of matrices for large problems, as described further below. In (28) the parameter *υ* controls the relative fit to data and dynamics. If both the data and dynamical error covariance are correct, *υ* should be one. We introduce this extra parameter in OTIS to allow for the possibility that covariances are misspecified. Once all of the matrices in (28) have been computed, trial solutions with different values of *υ* are readily computed, so this parameter can be tuned to adjust the relative fit to data and dynamics.

Clearly we would like to choose the subset of *N* locations for the calculated representers so that we do not loose important details in the approximate inverse solution (25). The best way to ensure this is far from clear, and some experimentation with different subsets of calculated representers may be necessary. We consider an example of this process in section 5.

### b. Processing of altimeter data

Here we describe processing steps used to prepare TOPEX/Poseidon altimeter data for the inversion. The starting point is a version of the PATHFINDER database (Koblinsky et al. 1999) with no tidal corrections applied (B. Beckely 1999, personal communication). To prepare the altimeter data for inversion, we first low-pass filter the altimeter data along each ground track, and decimate from the original spacing of ∼7 km (measured by the altimeter at a 1-s interval) to a spacing comparable to the grid resolution. Simple averaging along track is reasonable since all data in the window are observed at essentially the same time (compared to the tidal cycle) and at the same location (compared to the grid resolution). During the averaging process a simple median filter is used to remove outliers.

We allow for *N*_{c} constituents and use the superscript *l* to denote harmonic constants for individual constituents. The altimeter data at time *t* and location **x** can then be written as

where *ζ*^{l}_{SE}(**x**) is the solid earth tide (e.g., Hendershott 1977); ɛ represents error (for purposes of estimating the tides this includes all nontidal oceanography as well as measurement noise), and *α*^{l}(*t*) = (1 + *f*^{l}(*t*)) exp{*i* [*ω*^{l}(*t* − *t*_{0}) + *V*^{l}(*t*_{0}) + *p*^{l}(*t*)]} gives the slowly modulated periodic variation for constituent *l.* In this expression *ω*^{l} is the tidal frequency, *V*^{l}(*t*_{0}) is the Greenwich phase, and *f*^{l}, *p*^{l} are nodal corrections for amplitude and phase.

We compute *ζ*^{l}_{SE} from the TPXO.3 model of Egbert (1997), using methods described in Ray and Sanchez (1989). As with the loading and self-attraction *ζ*^{l}_{SAL} needed for forcing the forward model, *ζ*^{l}_{SE} is essentially a smoothed version of the ocean tide *ζ*^{l}. A similar (and reasonably accurate) solid earth tide estimate results from any of the recent global ocean tide models considered by Shum et al. (1997). We use the estimated *ζ*^{l}_{SE} to correct the altimeter data for solid earth tides, and at the same time we subtract the prior model *ζ*_{0}:

Finally, averages (over all times at each fixed location) are subtracted from the data to eliminate geoid errors. At a fixed data location **x**_{k} the result is a time series sampled at times *t*_{n}, *n* = 1, … , *T*_{k}, which we denote by **d**_{k}. This *T*_{k} dimensional vector satisfies

where **z**_{k} is the *N*_{c} dimensional complex vector of corrections to the prior solution at **x**_{k} with elements *z*^{l}_{k} = *ζ*^{l}(**x**_{k}) − *ζ*^{l}_{0}(**x**_{k}), *l* = 1, … , *N*_{c}. Here 𝗔_{k} is the *T*_{k} × *N*_{c} matrix with *l, **n* element

For TOPEX/Poseidon, time series at most locations are now long enough to separate all major constituents. It is thus possible to estimate the vectors of harmonic constants **z**_{k} by least squares fitting of (31) (see appendix C). The resulting along-track estimates of *z*^{l}_{k} can then be used as data for constituent *l* at location **x**_{k} in the reduced basis inversion described above.

### c. Time domain data

For short time series harmonic analysis of data from a single location into individual constituents may not be possible. As an extreme example, consider the case of shipborne acoustic Doppler current profiler measurements of currents, for which there will generally be only a single observation at any fixed location. Although any sort of harmonic analysis would obviously be impossible, such data could still be assimilated into a multiple constituent tidal model. In effect, the spatial structure imposed by the dynamics allows data from multiple locations but different times to be combined in a rational way, so that individual constituents can still be separated (at least approximately) over the full model domain.

The major complication with time-domain data is that the equations for the representer coefficients for all *N*_{c} constituents are coupled. Matrices required for the inversion are thus much larger (by a factor of *N*_{c} for both rows and columns), making computation of the representer coefficients for large datasets a challenge. Some additional steps and approximations are thus required to keep computations tractable. A second minor complication with the time domain inversion is that the data are intrinsically real, while the tidal constituents (and the dynamical equations) are most simply expressed in the complex domain. This detail is easily dealt with, though it complicates notation, and leads to a further increase in storage requirements (by a factor of 2) for some matrices. The time-domain tidal data inversion problem is treated in some detail in EBF, along with the more general case of dynamical errors that are correlated between constituents of a fixed species. This correlation also leads to coupled equations for the representer coefficients, even for harmonically analyzed data. Here we only summarize the key points for the simpler case of time domain data with no interconstituent correlation of dynamical errors.

As we show in appendix C the data vector for a single site can be reduced to at most *M*_{k} = min(2*N*_{c}, *T*_{k}) elements. This data reduction is accomplished by a rotation in the data space with a 𝗤𝗥 decomposition, and is thus unconditionally stable (and hence safe even if constituents are poorly separated). The *M*_{k} dimensional reduced data vector **d̃**_{k} still satisfies an equation of the form given in (31), with the complex matrix 𝗔_{k} replaced by the *M*_{k} × *N*_{c} complex matrix 𝗕_{k} defined in appendix C. Let 𝗥^{l}, 𝗣^{l}, and 𝗘^{l}, *l* = 1, *N*_{c} be as for the single constituent case, and define the block diagonal matrices 𝗣 = diag(𝗣^{l}, *l* = 1, *N*_{c}) and 𝗘 = diag(𝗘^{l}, *l* = 1, *N*_{c}). Let ** β** be the representer coefficients for all constituents, and let

**d̃**be the reduced data for all sites:

In appendix C we construct a sparse *M* × *M* matrix 𝗕 (where *M* = Σ^{K}_{k=1 }*M*_{k}) from the elements of the matrices 𝗕_{k}, *k* = 1, *K* such that the full reduced data vector satisfies:

The real matrix **Σ**^{−1}_{e}𝗙 is analogous to the complex matrix **Σ**^{−1}_{e}𝗣𝗘 for the one constituent case. The solution coefficient vector ** β** can be found in a similar fashion, by computing the singular value decomposition (SVD) of

**Σ**

^{−1}

_{e}𝗙 = 𝗪𝗦𝗤

^{T}, solving for the real coefficients

**c**

converting **c** to a complex vector **c** = **c**_{R} + *i***c**_{I}, and then computing the full vector of representer coefficients (for all constituents) as ** β** = 𝗘

**c**.

Note that the matrix 𝗙 is *M* × (2*N*_{c}*N*), and for large datasets (and many tidal constituents) this matrix can easily be too large to store in RAM or to SVD all at once. Two modifications to the basic procedure are used in OTIS to make multiconstituent inversion of large datasets feasible. First, we reduce the sizes of the matrices by truncating the matrices 𝗘^{l} to include only the largest *N*_{T} < *N* eigenvalues of 𝗥^{l}. This reduces the size of 𝗙 to *M* × (2*N*_{c}*N*_{T}). Second, we compute the SVD with a two stage procedure, which allows an efficient out of core algorithm. Further details are given in appendix C.

### d. Posterior error calculation

Posterior errors for the inverse solution can be estimated using the Monte Carlo approach described in Dushaw et al. (1997). Although these error bars will only be as reasonable as the prior errors of both the data and dynamics, they can provide at least some guidance in assessing solution reliability. See Dushaw et al. (1997), Egbert and Ray (2000), and Ray et al. (2001) for example uses of these error bars. The calculation is organized in four steps:

- Generate a series of
*I*random forcings and boundary conditions with the assumed dynamical error covariance, and solve the forward problem for this forcing. That is, for realization*i*we compute:where(36)*δ***u**_{i}^{−1}**Σ**^{1/2}_{f}**w**_{i}**w**_{i}is a vector of unit variance uncorrelated complex Gaussian pseudorandom numbers. In the representer calculation, smoothing with covariance**Σ**_{f}is accomplished by pseudo–time stepping a diffusion equation a total of*T*steps (EBF). Multiplication by**Σ**^{1/2}_{f}is thus accomplished by stopping the smoothing after*T*/2 steps. It is readily verified that the resulting smoothed random forcing error fields**Σ**^{1/2}_{f}**w**_{i}have covariance**Σ**_{f}. We estimate an adjustment factor for the assumed dynamical error**Σ**_{f}, which makes the overall (spatially averaged) variance of the synthetic realizations consistent with the magnitude of the model change for the actual inverse solution. The adjustment factor is computed using the difference between the prior solution and the inverse solution. If the dynamical error magnitude assumptions are correct, then the adjustment factor should be near one. For each realization evaluate the tidal error field

*δ***u**at the data locations**x**_{k}and add measurement noise (scaled with the assumed data error) to compute a synthetic data vector**d**_{i}. Again, an adjustment factor is computed to make the synthetic data error consistent with the actual inverse solution residuals. If data error assumptions are correct, then the adjustment factor should be near one. The synthetic data vector then substituted into (28) or (35) to solve for the representer coefficients*β*_{i},*i*= 1,*I.*For each

*β*_{i}, use (21)–(23) to compute the inverse solutions**û**_{i}=**u**_{0}+*δ***u**_{i}(including elevations and transports).- Using the actual synthetic solutions
**u**_{i}=**u**_{0}+*δ***u**_{i}and the inversion results**û**_{i}, posterior error variances and covariances for any functional of the tidal fields can then be calculated. For example the pointwise variance of the estimated elevation for constituent*l*is computed as Note that this procedure can also be applied to calculation of error bars for more general (possibly nonlinear) functionals [e.g., spherical harmonic coefficients of the tides (Ray et al. 2001) or integrated energy flux into shallow seas (Egbert and Ray 2000)].

## 5. Example of T/P data inversion for tides in the Indonesian Seas

Both the diurnal and semidiurnal tides in the semienclosed Indonesian Seas are spatially complex and of large amplitude. Hatayama et al. (1996) review tidal studies and present results from a new purely hydrodynamic model of the area. Tidal currents are large enough to be of great practical importance to navigation in many locations, and there is evidence that the tides may play a crucial role in interchange processes between the Pacific and Indian Oceans (Hatayama et al. 1996). Because the tides are greatly influenced by the complex topography of the area they are quite difficult to model accurately with a purely dynamical model (Mazzega and Berge 1994).

TOPEX/Poseidon altimeter data provide useful constraints on tidal elevations, at least in the larger seas within the Indonesian archipelago. Elevation maps for this area are included in most of the T/P-based global tidal solutions discussed by Shum et al. (1997). However, there are great differences between all of the various solutions, and as we shall show below for the TPXO.3 solution of Egbert (1997), misfits to estimates of tidal harmonic constants computed along the altimeter track can be quite large, frequently exceeding several tens of centimeters. For these global ocean tide models, grid resolution is too coarse to expect high accuracy, given the complex and often shallow bathymetry of the Indonesian seas. Higher-resolution assimilation or inverse modeling is called for in this topographically complicated area.

Mazzega and Berge (1994) assimilated T/P data along with tide gauges to model tidal fields in the Asian semienclosed seas, particularly in the South China Sea and in the Indonesian Seas (10°S–15°N, 105°–135°E). Though elevations in this model were shown to be a significant improvement compared to Schwiderski (1980), the relatively coarse ½° resolution used was insufficient for estimating details of the tides in this area (Hatayama et al. 1996). The coarse resolution and limited-area used by Mazzega and Berge (1994) were dictated by the available computer resources. Despite continued rapid progress in computer technologies, hydrodynamic modeling and assimilation in such a complex area remains a challenging problem, especially given the growing volume of data available for assimilation. The Indonesian Seas thus provide a good illustration of the effectiveness of the efficient GI approach developed in this paper. Here we concentrate on demonstrating the advantages and effectiveness of the approach, rather than providing a detailed interpretation of the dynamics or consequences revealed by the inversion results.

We consider a model domain with longitude and latitude limits 21°S–15.6667°N, 95°–165°E (70° × 36.6667°), shown in Fig. 2. The grid used is 420 × 240, with ⅙° resolution. Bathymetry was computed by averaging the Gtopo30 topography database (Smith and Sandwell 1997) onto the *ζ* nodes of the C grid. Eight principal tidal constituents (*M*_{2}, *S*_{2}, *K*_{1}, *O*_{1}, *N*_{2}, *P*_{1}, *K*_{2}, *Q*_{1}) are modeled.

The inversion scheme described in the previous sections was applied to 232 cycles of T/P altimeter data from a modified version of the pathfinder database (Koblinsky et al. 1999) with no tidal corrections applied (B. Beckeley 1999, personal communication). For this example we used altimeter data for 5880 data sites, including all crossover points in the model domain, and each fifth location between crossover points, with a spacing of ∼28 km. The data locations are illustrated in Fig. 3. A prior solution was obtained by time stepping the nonlinear barotropic shallow water equations (including advection, nonlinearities in the continuity equation, a quadratic parameterization of bottom friction, and lateral eddy viscosity). Sea surface elevations from TPXO.3 (Egbert 1997) were used as the open boundary conditions. As described in section 4, the T/P data were then low-pass filtered along track before sampling at the data sites, and equilibrium long-period tides, a correction for radial deformation, the prior solution, and the time series average were subtracted from the altimeter data at each data location. Finally, data were harmonically analyzed for the eight modeled constituents, as described in appendix C.

Using Eq. (20) 312 representers at the locations shown by the filled symbols in Fig. 3 were calculated. The dynamical error covariance **Σ**_{f} was defined following the considerations outlined in EBF, using the prior solution to estimate the spatially varying magnitudes of errors in the momentum equations. The correlation length scale for the dynamical errors was set to 200 km (approximately eight grid cells). The continuity equation was assumed to be exact. An example of an elevation data representer for an observation of *M*_{2} in the South China Sea is given in Fig. 4.

The bulk of the computation required for implementation of the GI approach lies in the representer calculation. Computational resources required to calculate representers for all eight constituents are compared in Table 1, for several different computer platforms, and for the time stepping and matrix factorization approaches. Factoring the coefficient matrix for the wave equation dramatically reduces the CPU time required for the representer calculation.

The direct factorization approach allows for a decrease in computational time (relative to the time-stepping approach used for the representer calculation of EBF) by a factor of approximately 3000. This allows solution of moderate size problems on a standard workstation that would not be practical with the time-stepping approach, even on a small supercomputer like the CM-5E. Note, however, that the required RAM increases rapidly with grid size. Storage of the banded matrix **A** requires (3*m* + 4)*N*_{oc} complex words, where *N*_{oc} is the number of ocean nodes (excluding land). For our 420 × 240 grid, 446 Mb of RAM are required for single precision computations. For double precision calculations close to 1 Gb of RAM would be required. This is near the limit of RAM for a desktop workstation at this time. Computations on significantly larger grids would require use of a different class of computer, or modifications to the algorithms to allow for blocked out-of-core factorization schemes.

To compute the inverse solution, we used the single constituent reduced basis approach [see Eqs. (25)–(28)] to calculate the representer coefficients ** β**. The inverse solution was then calculated using Eqs. (21)–(23). Figure 5 shows elevations for the prior and inverse solution for the

*M*

_{2}and

*K*

_{1}tidal constituents. The most significant changes from the prior solution are in the shallow South China, Java, and Arafura Seas, in the Gulfs of Carpentaria and Thailand, and in the Sulu Sea. The amplitude of change in the elevation fields is as great as 60 cm for

*M*

_{2}, and 40 cm for

*K*

_{1}.

Figures 6 and 7 show in-phase and quadrature tidal currents for *M*_{2} and *K*_{1}, for the inverse and selected parts of the prior solutions. We show only the central part of the model domain, where seas are shallowest and currents largest. Assimilating T/P data results in significant changes in *M*_{2} tidal currents in the Java, Timor, and Arafura Seas, and along the northwest coast of Australia, where currents in the inverse solution reach a few meters per second. (Note that plotted vector magnitudes are limited to 1 m for clarity.) Fitting the altimeter data also results in significant phase shifts of currents in the Timor Sea, and of flow into the Sulu Sea and through Makassar Straight, one of main connections between the Pacific and Indian Oceans (Hatayama et al. 1996).

Major changes in the *K*_{1} diurnal currents occur in the adjacent parts of the South China and Java Seas, in the Arafura Sea, and in the Gulf of Carpentaria, and also in the Molucca, Halmahera, and Seram Seas (to the northwest of New Guinea).

In Fig. 8 we plot the difference between the inverse solution elevation and estimates of *M*_{2} harmonic constants computed along the T/P ground track. For comparison we also plot residuals computed using the TPXO.3 global solution. For the global solution, residuals are much larger, and spatially coherent (i.e., nearby residuals all tend to have real and imaginary parts of the same sign). For the inverse solution residuals mostly look like noise, with residuals on crossing tracks generally having opposite signs. This implies that the altimeter data is not consistent with a dynamically realistic tidal solution, so that a significantly better fit is not possible.

Table 2 summarizes the root-mean-square residual magnitude for the four main tidal constituents as a function of water depth for the prior and inverse solutions. Results are also given for TPXO.3 and for two other more recent global solutions—GOT99 (Ray 1999) and CSR4.0 (R. Eanes 1999, personal communication). The inversion procedure dramatically improves fit to the data compared to the purely hydrodynamic prior model in all depth ranges, and for all constituents. Misfits for all solutions increase in shallower water. In the shallowest depth range (*H* < 500 m) the coarse (½° resolution) global models generally fit the data rather poorly, with rms misfit for *M*_{2} of the order of 10 cm, although GOT99 fits the altimeter data significantly better in these shallow depths. The higher resolution inverse solution provides by far the best fit in shallower water. In deeper water our regional inverse solution generally does not fit the altimeter data significantly better than the global solutions (especially compared to GOT99).

Although the deep water open boundary conditions obtained from the global inverse solution TPXO.3 are reasonably accurate, there are still significant discrepancies between the prior solution and the altimetry data in many shallow areas. This demonstrates that there are significant errors in the dynamical equations used to model the tide in this area. The two largest sources of error are almost certainly in the bathymetry and in the representation of dissipation. Discrepancies in this area between two widely available digital global topographic databases [ETOPO5 (National Geophysical Data Center 1992) and Gtopo30 (Smith and Sandwell 1997)] reach ∼100 m over large areas in the shallows. It is not clear which database is most accurate; most likely both contain significant errors in this topographically (and politically) complex area. For the prior model we assumed the quadratic parameterization for bottom friction of (11). This is at best an approximation, which must limit solution accuracy in shallow seas where dissipative terms are significant. Note also that there are areas of steep topography within our domain, where generation of internal tides by barotropic tidal flow may be significant (e.g., Sjoberg and Stiegbrandt 1992; Egbert and Ray 2001). These baroclinic processes, which are not accounted for in the shallow water equations or our parameterization of dissipative processes, can extract significant energy from the barotropic tide.

In the last column of Table 2 we give the average least squares estimation errors for the along-track harmonic constants. These give some indication of noise levels in the data. In shallow water misfits to the inverse solution are slightly larger than the harmonic constant standard errors, suggesting slightly better fits to the data might be possible. However, some small-scale features in the along-track tidal elevation harmonic constants may be associated with baroclinic tides that are phase locked to the barotropic tide (Ray and Mitchum 1997), and this part of the signal should not be fit by a barotropic model. In deeper water fit of the inverse solution is about as good as possible. Overall, we conclude that the inverse solution fits the altimeter data significantly better than previous solutions, and probably about as well as can be justified, given the level of noise in the data.

With the reduced basis approach, it is important to choose the locations for the calculated representers so that we do not loose important details in the approximate inverse solution. For the Indonesian Seas solution, we experimented with five different subsets of representers, which are indicated by the different symbols in Fig. 3. The first subset consists of 90 representers spread uniformly on a subset of crossover points, denoted with the filled squares. The second set includes an additional 90 representers marked by filled circles (the remaining crossover points). The remaining subsets have an additional 60, 30, and 42 representers, which are denoted by diamonds, triangles, and inverted triangles, respectively. These are placed between crossover points in shallow seas and areas of complex bottom topography. The final sets of representers were chosen after examination of residual plots (such as shown in Fig. 8) for preliminary versions of the inverse solution. For the results discussed above we used the final full set of 312. Figures 9a,b shows along-track rms misfit versus the number of representers *N* used for the solution, for the two shallower depth ranges (*H* < 500 m, 500 m < *H* < 2500 m). The value at *N* = 0 corresponds to the prior model misfit. For deeper areas misfits are reduced to near the lowest levels with approximately 100 representers. For the shallowest depth, range misfits are still falling with the addition of more representers, suggesting that small additional improvements may be obtained, particularly if the along-track misfit plot of Fig. 8 is used for guidance in choosing representer locations. Alternatively, the conjugate gradient approach described in Egbert and Bennett (1996) can be used to refine the solution. However, as discussed above, given possible complications due to surface expression of the internal tides, it is not clear that better fits are actually justified.

Figures 10a,b show posterior errors for the *M*_{2} and *K*_{1} elevations estimated using the Monte Carlo simulation approach and Eq. (37). In the open ocean, errors are less than 1 cm. Larger errors (4 cm) occur in the shallower water, especially in the coastal areas where tidal amplitudes are greatest, and residuals in Fig. 8 are largest. Similar pattern of error amplitudes are observed for both diurnal and semidiurnal constituents.

## 6. Summary and conclusions

In principal the GI approach has a number of significant advantages over simpler approaches to tidal data assimilation such as nudging (Kantha 1995; Matsumoto et al. 1995), or optimal interpolation (Mazzega and Berge 1994). The GI approach allows complete control over conditioning of the inverse solution, provides a natural framework for rational specification of errors in the dynamical equations and data, and is optimal for the assumed covariances. Also, GI allows for the computation of posterior error bars, and testing of hypotheses (e.g., Bennett 1992; Egbert 1997). However, the rather heavy computations required for GI are a significant enough disadvantage to have limited the use of this approach in practice. Here we have developed and implemented an efficient linearized scheme, which makes routine application of GI to global and regional barotropic tidal modeling quite feasible.

To make the inversion efficient, we factor the coefficient matrix for the frequency domain elevation wave equation, derived from the linearized shallow water equations by eliminating velocities. Once this matrix has been factored, the repeated solution of the linearized tidal equations required to calculate representers is very fast. This direct factorization approach is more than a factor of 1000 faster than the time step representer calculation used by EBF to calculate representers for a global tidal inverse solution. Combined with a reduced basis representer approach this solution scheme makes application of the GI approach to regional-scale tidal modeling problems (with approximately 10^{5} grid nodes and 10^{4} data) practical even on a desktop workstation.

To keep memory requirements for storing the factored coefficient matrix reasonable, it is necessary to eliminate velocities from the coupled first-order system of shallow water equations (12)–(13) to derive a wave equation in the scalar elevation field. A key feature of our approach is that we keep the penalty functional defined in terms of the first-order equations, which embody conservation of momentum and mass of the ocean fluid, so we can allow for errors in these original equations. If the penalty functional were instead formulated directly in terms of the wave equation, velocities could not be computed from the inverse elevation solution without assuming that the momentum equations (12) are exact. In effect all of the error would then be assumed to be in Eq. (13). In fact, this equation is just a simple statement of mass conservation, which (with the equations expressed in terms of volume transports) involves no unknown or uncertain parameters. With a reasonable interpretation of the transports and elevations on the C grid as averages over cell sides and the cell area, respectively, there is not even any grid truncation error in the numerical approximation of this equation. In estimating currents it would be unreasonable to not enforce mass conservation, but to require an exact fit to the momentum equations, which involve uncertain bathymetry, and approximate parameterizations of dissipation and tidal loading and self-attraction. The proper choice of covariance has important implications for some applications. Egbert and Ray (2001) show that empirical estimates of energy dissipation obtained without enforcing mass conservation are too noisy to be interpreted, while stable and informative estimates result when fit to (13) is enforced.

With the approach developed here it is possible to enforce mass conservation exactly, leaving all error in the momentum equations. It is also possible to assimilate tidal currents, even though all dynamical calculations are done strictly in terms of elevations. This general approach could be usefully applied to other assimilation problems where some model variables can be eliminated to simplify dynamical calculations, and then calculated prognostically from those variables retained. The basic idea is embodied in Eq. (18): by writing explicit expressions for the prognostic calculation equations, together these with the reduced dynamical equations, it is straightforward to derive expressions for the appropriate adjoints [i.e., Eq. (19)]. These show how to compute the forcing error for the original coupled system, using adjoints of the reduced system and of the prognostic calculations.

Our efficient representer calculation scheme has been combined with programs for generating grids, boundary conditions, tidal forcing, dynamical error covariances, and altimetry datasets into a relocatable system for inverse barotropic tidal modeling. We have illustrated the functionality of the system by constructing a new regional-scale inverse solution for the main eight tidal constituents in the complex Indonesian Seas. The solution fits the altimetry data significantly better than previously proposed coarser-scale global tidal models. Our results suggest that with careful implementation assimilation based on a rigorous GI approach can be quite efficient.

## Acknowledgments

We thank Richard Ray for providing tidal loading and self-attraction corrections, and for discussions, and Brian Beckely for the no-tides pathfinder database. This work was supported by the Office of Naval Research through Grant N00014-9410926 and by the National Science Foundation through Grant OCE-9633527.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{2}ocean tide.

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Numerical Implementation Details

We consider operators 𝗚 (gradient), 𝗗 (divergence), and 𝗖 (multiplication by **Ω**^{−1}) as defined in section 3, Eq. (17). The nodes *ζ*_{ij}, *U*_{ij}, and *V*_{ij} are defined as in Fig. 1, section 3. Depth *H*_{ij} is defined on the *ζ* nodes. Let *M*^{U}, *M*^{V}, and *M*^{ζ} be the water/land masks (0 for land nodes, 1 for ocean nodes). The node *U*_{ij} or *V*_{ij} is land if any adjacent *ζ* node is. In particular masks are zero for coastal boundary nodes. Then operators 𝗚, 𝗗, and 𝗖 can be expressed explicitly on the C grid as

Here *a* is the radius of the earth, *g* is gravitational acceleration, (*θ, **ϕ*) are the latitude and longitude. The 2 × 2 matrices **Ω**^{U}_{ij} and **Ω**^{V}_{ij} are defined by (14) at the *U*_{ij} and *V*_{ij} nodes, and the following components of the inverse matrices are used in (A3):

The operator 𝗔: *ζ* → *ζ* from (17) is defined as a concatenation of 𝗗, 𝗖, and 𝗚, plus the obvious scalar multiplication by *iω.*

Land nodes are excluded from transforms (A1), (A2), and (A3) by the masks. In particular, the gradient operator is zero on all boundary *U* and *V* nodes. The operators 𝗔 and 𝗖 require special definitions on boundary nodes. Here 𝗖 is identity on *U*/*V* boundary nodes:

and 𝗔 is identity on *ζ* boundary nodes:

It is easily verified from (18) that with these conventions, the solution on any boundary node (say *U*_{ij}) will satisfy *U*_{ij} = (**f**_{U})_{ij}, and similarly for the *V* and *ζ* nodes. With these conventions general inhomogeneous boundary conditions at the coast and open boundaries can thus be enforced. Note that this general form for the boundary conditions is required because we allow for errors in all boundary conditions (in particular in the no-flow boundary condition at the coast). These conventions define a consistent treatment of boundary conditions both for the wave equation, and for the computation of volume transport. Note also that the special treatment of boundary nodes also defines [through the discrete adjoint (19)] the correct adjoint boundary conditions for the representer calculation.

### APPENDIX B

#### Global Grids

For a global grid that is cyclic in longitude the coefficient matrix for the wave equation is not exactly banded anymore. To allow for this case, we partition the coefficient matrix into four submatrices, where the largest is banded. In general for any matrix 𝗔 and vectors **b**, **x**, partitioned as

the solution to 𝗔**x** = **b** satisfies

Let *n* be the number of longitude divisions, or columns in the grid, and *m* be the number of latitude divisions, or rows in the grid. First we can partition the coefficient matrix for the wave equation 𝗔 into *n* rows and *n* columns as in (B1), where

The largest matrix 𝗔_{11} corresponds to the full coefficient matrix 𝗔, except all columns and rows involving the rightmost column of the grid. Here 𝗔_{11} has the usual banded structure and can be factored as in section 3. The matrix 𝗔_{12} consists of *m* columns *α*_{k}, *k* = 1, *m.*

Thus to solve the wave equation for the global case, we do the following.

(a) Solve 𝗔_{11}**y** = **b**_{1} for **y** and 𝗔_{11}**z**_{k} = *α*_{k} for **z**_{k}, *k* = 1, *m* using the factorization of 𝗔_{11} (the factorization need be done only once). Collect **z**_{k}, *k* = 1, *m* into the matrix 𝗭 = |**z**_{1}, **z**_{2}, … , **z**_{m}|.

(b) Compute **w** = **b**_{2} − 𝗔_{21}**y** and the square *m* × *m* matrix 𝗘 = 𝗔_{22} − 𝗔_{21}𝗭.

(c) Solve the full *m* × *m* system 𝗘**x**_{2} = **w**.

(d) Compute **u** = **b** − 𝗔_{12}**x**_{2}.

(e) Solve 𝗔_{11}**x**_{1} = **u**, again using the factorization of 𝗔_{11}.

### APPENDIX C

#### Linear Algebra and Data Reduction

We first consider the initial reduction of the time domain data. We rewrite (31) explicitly in terms of real and imaginary parts

where the overbars are used to denote real versions of complex matrices and vectors. Assuming *T*_{k} ≥ 2*N*_{c}, we next compute the 𝗤𝗥 decomposition of the real matrix 𝗔_{k} = 𝗤𝗕_{k}. Here 𝗤 is of size *T*_{k} × 2*N*_{c} with orthonormal columns and 𝗕_{k} is an upper-triangular real matrix of size 2*N*_{c} × 2*N*_{c}. Multiplying both sides of (C1) by 𝗤^{T} reduces the data vector to 2*N*_{c} components

In the last step we have converted the real matrix 𝗕_{k} back to complex via 𝗕_{k} = 𝗕_{k1} − *i*𝗕_{k2}, where 𝗕_{k} = [𝗕_{k1 }𝗕_{k2}] is partitioned into two real submatrices, each with *N*_{c} columns. With the transformation of (C2) we have reduced the data vector to 2*N*_{c} components, discarding the part of the data that cannot be fit by any *N*_{c} constituent tidal model. Using the usual least squares (LS) theory, we estimate the data error variance from the mean-square residual magnitude *σ̂*^{2}_{k} = ‖**d**_{k} − **d̃**_{k}‖^{2}/(*T*_{k} − 2*N*_{c}). Assuming that the original data errors **ɛ**_{k} are uncorrelated with variance *σ*^{2}_{k}, the transformed error vector ɛ_{k} will have the same simple covariance, since multiplication by 𝗤^{T} corresponds to a rotation in the data space. We make this simplifying assumption in OTIS, and use *σ̂*^{2}_{k} to estimate the error variance for the 2*N*_{c} elements of **d̃**_{k}. For each data location **x**_{k} this initial data reduction procedure yields a reduced data vector **d̃**_{k}, temporal sampling structure matrix 𝗕_{k}, and data error variance *σ̂*^{2}_{k}. These can then be used along with the calculated representers to solve simultaneously for all *N*_{c} constituents as outlined in section 4. Note that if the time series at **x**_{j} is shorter than 2*N*_{c}, there is no need for this initial reduction of the data vector, and the matrix 𝗔_{k} can be used in place of 𝗕_{k}. However in this case an independent estimate of the data error variance *σ*^{2}_{k} is required.

For long enough time series constituents can be separated at each location, allowing the more efficient single constituent inversion approach to be used. To complete harmonic analysis of the time series **d**_{k}, we can simply solve the triangular system **d̃**_{k} = 𝗕**z**_{k}. The result will be the usual LS estimate of the real and imaginary parts of the *N*_{c} tidal harmonic constants. Error variances for the elements of **z**_{k} are obtained in the usual way for the LS estimates from the diagonal elements of *σ̂*^{2}_{k}(𝗕^{T}_{k}𝗕_{k})^{−1}. These estimates of error variances for each constituent are then used to define the diagonal data error covariance **Σ**^{2}_{e}.

Next, we consider further details on matrix manipulations required for the multiple constituent inversion. Let **B**_{k} = |**b**^{1}_{k }**b**^{2}_{k} · · · **b**^{Nc}_{k}| denote the *N*_{c} columns of the complex matrices defined above. Define the large sparse blocked matrix

With 𝗕 so defined, it is then straightforward to check that (34) holds.

Finally we consider computation of the SVD of the potentially very large matrix **Σ**^{−1}_{e}𝗙, where 𝗙 is defined in (34). First note we can reduce the number of columns in each of the matrices 𝗘^{l} by keeping only the eigenvectors of 𝗥^{l} associated with the largest *N*_{T} eigenvalues. Typically *N*_{T} can be approximately half of *K* with little effect on the final approximate inverse solution. With this reduction **Σ**^{−1}_{e}𝗙 is *M* × 2*N*_{c}*N*_{T}. Since *M* can be as large as 10^{5}, this matrix can still be too large to allow direct computation of the SVD in core. For OTIS we use a two-stage procedure. In the first stage we compute a 𝗤𝗥 decomposition **Σ**^{−1}_{e}𝗙 = 𝗪_{F}𝗤_{F}, where 𝗤_{F} is an upper-triangular matrix of size 2*N*_{c}*N*_{T} × 2*N*_{c}*N*_{T}. Computation of 𝗤_{F} and 𝗪^{T}_{F}**d̃** can be accomplished efficiently with a blocked out-of-core algorithm (e.g., Golub and Van Loan 1989). Note that this scheme does not directly compute or store the large orthogonal matrix 𝗪_{F}. In the second stage we compute the SVD of 𝗤_{F} = 𝗪_{2}𝗦𝗤^{T}. Then **Σ**^{−1}_{e}𝗙 = (𝗪_{F}𝗪_{2})𝗦𝗤^{T} = 𝗪𝗦𝗤^{T} is the SVD of the large matrix **Σ**^{−1}_{e}𝗙. Note that the product 𝗪^{T}**d̃** = 𝗪^{T}_{2}𝗪^{T}_{F}**d̃** [required for computation of the representer coefficients using (35)] is easily computed without actually storing 𝗪. It is useful to compute and store this large orthogonal matrix for the posterior error calculation only. This last matrix multiplication to compute 𝗪 can also be easily accomplished with an out-of-core algorithm.

## Footnotes

*Corresponding author address:* Dr. Gary D. Egbert, College of Oceanic and Atmospheric Sciences, Ocean Admin. Building 104, Oregon State University, Corvallis, OR 97331-5503. Email: egbert@oce.orst.edu