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 105 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 105 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

 
uf0
(1)

and with a K-dimensional vector d of tidal data

 
du
(2)

In (2) 𝗟 = [L1, … , LK] 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

 
formula

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

 
formula

where u0 = 𝗦−1f0 is the exact solution of (1) and the functions rk are the representers (Yosida 1980) of the data functionals defined by Lk, k = 1, K. Representers can be calculated by first solving the adjoint of the dynamical equation

 
αkk
(5)

where Δk is the averaging kernel for the data functional Lk (an impulse at xk for a point observation at this location), and then solving the forward equation

 
rkΣfαk
(6)

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

 
Σeβdu0
(7)

where 𝗥 is the representer matrix with elements

 
RjkLjrk
(8)

In the most straightforward application of the representer approach, (5) and (6) are solved K times each for the functions rk, 𝗥 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

 
formula

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 f0.

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 = f0 + 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)

 
FcDvHU
(11)

where v is the total velocity vector (in particular, including all tidal constituents), and the value of the nondimensional parameter cD 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 u0 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 κ = (cDv‖/H). Several refinements are possible. The time domain equations can be linearized around the solution u0, 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., M2, S2, K1, O1) 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 xk. In this case the forcing for the representer calculation (8) (i.e., Lk) corresponds to an impulsive forcing at xk. 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

 
formula

where

 
formula

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

 
UgHΩ−1ζΩ−1fU
(15)

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

 
gHΩ−1ζiωζΩ−1fUfζ
(16)

Solution of (1) can thus be accomplished by solving (16) for ζ, then using the result in (15) to directly calculate U. This provides our basic scheme for solving the shallow water equations.

a. Numerical implementation

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

 
formula

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 ), maps from ζ nodes to ζ nodes and represents the second-order operator of the frequency domain wave equation.

Fig. 1.

Grid conventions used in OTIS. An Arakawa C grid is used for all dynamical calculations. Volume transports U and V are specified on grid cell edges, and are interpreted to be the average volume transport over the cell edge. Elevations are interpreted as the average over the cell, and are given at the center. Boundary conditions at the coast are specified on the U and V nodes. Open boundary conditions are given by specifying the elevation ζ for open boundary edge cells, or transports on edge U or V nodes

Fig. 1.

Grid conventions used in OTIS. An Arakawa C grid is used for all dynamical calculations. Volume transports U and V are specified on grid cell edges, and are interpreted to be the average volume transport over the cell edge. Elevations are interpreted as the average over the cell, and are given at the center. Boundary conditions at the coast are specified on the U and V nodes. Open boundary conditions are given by specifying the elevation ζ for open boundary edge cells, or transports on edge U or V nodes

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

 
formula

Arbitrary inhomogeneous boundary conditions at the coast and open boundaries are readily implemented through appropriate definitio of 𝗚,𝗗,𝗖,𝗔, fζ, and fU 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:

 
formula

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 (3m + 4)nm 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 xk, k = 1, K. In this case the averaging kernel for data functional Lk takes the form of an impulse in the ζ component at xk. 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 UV nodes as in (18), this implies

 
formula

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

 
rζ−1Σf−1*k
(20)

where Δk is a unit magnitude impulsive forcing at xk. The transport component of the representer rU 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:

  1. Using the LU factorization of 𝗔, solve the conjugate transpose system with impulsive forcing at xk: 𝗔*αζ = Δk. The solution is defined on the ζ nodes.

  2. Apply the local operators 𝗗* and then 𝗖*. These are the transpose of 𝗗 and 𝗖, respectively, and map αζ to αU, defined on the U/V nodes.

  3. Apply the covariance smoother Σf to αU.

  4. 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.

  5. 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., Rjk is computed by evaluating rk at xj) 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

 
formula

This is exactly the forcing error required to bring the inverse solution û = u0 + δ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

 
δζ−1δf
(22)

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

 
δUδζδf
(23)

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 xk 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

 
formula

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 105, 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 Pkj = Lkrj, 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)

 
formula

For u defined by (25) we clearly have

 
formula

Also, EBF show

 
𝗦(uu0)*Σ−1f𝗦(uu0) = β*𝗥β

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

 
JβddβΣ−1edβββ
(26)

where d′ = d − 𝗟u0. 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

 
Σ−1e
(27)

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 si. It is simple to show (e.g., Parker 1994) that the minimizer of (26) can then be written as

 
β2υ−1d
(28)

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.

Once the representer coefficients β are computed these are used in (25) or (21)–(23) to compute the correction δu to the prior solution.

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 Nc 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

 
formula

where ζlSE(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 + fl(t)) exp{i [ωl(tt0) + Vl(t0) + pl(t)]} gives the slowly modulated periodic variation for constituent l. In this expression ωl is the tidal frequency, Vl(t0) is the Greenwich phase, and fl, pl are nodal corrections for amplitude and phase.

We compute ζlSE from the TPXO.3 model of Egbert (1997), using methods described in Ray and Sanchez (1989). As with the loading and self-attraction ζlSAL needed for forcing the forward model, ζlSE 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 ζlSE to correct the altimeter data for solid earth tides, and at the same time we subtract the prior model ζ0:

 
formula

Finally, averages (over all times at each fixed location) are subtracted from the data to eliminate geoid errors. At a fixed data location xk the result is a time series sampled at times tn, n = 1, … , Tk, which we denote by dk. This Tk dimensional vector satisfies

 
dkkzk
(31)

where zk is the Nc dimensional complex vector of corrections to the prior solution at xk with elements zlk = ζl(xk) − ζl0(xk), l = 1, … , Nc. Here 𝗔k is the Tk × Nc matrix with l, n element

 
formula

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 zk by least squares fitting of (31) (see appendix C). The resulting along-track estimates of zlk can then be used as data for constituent l at location xk 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 Nc constituents are coupled. Matrices required for the inversion are thus much larger (by a factor of Nc 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 Mk = min(2Nc, Tk) 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 Mk dimensional reduced data vector k still satisfies an equation of the form given in (31), with the complex matrix 𝗔k replaced by the Mk × Nc complex matrix 𝗕k defined in appendix C. Let 𝗥l, 𝗣l, and 𝗘l, l = 1, Nc be as for the single constituent case, and define the block diagonal matrices 𝗣 = diag(𝗣l, l = 1, Nc) and 𝗘 = diag(𝗘l, l = 1, Nc). Let β be the representer coefficients for all constituents, and let be the reduced data for all sites:

 
formula

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

 
formula

The real matrix Σ−1e𝗙 is analogous to the complex matrix Σ−1e𝗣𝗘 for the one constituent case. The solution coefficient vector β can be found in a similar fashion, by computing the singular value decomposition (SVD) of Σ−1e𝗙 = 𝗪𝗦𝗤T, solving for the real coefficients c

 
c
2υ−1T
(35)

converting c to a complex vector c = cR + icI, and then computing the full vector of representer coefficients (for all constituents) as β = 𝗘c.

Note that the matrix 𝗙 is M × (2NcN), 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 NT < N eigenvalues of 𝗥l. This reduces the size of 𝗙 to M × (2NcNT). 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:

  1. 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: 
    δui−1Σ1/2fwi
    (36)
    where wi 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/2f is thus accomplished by stopping the smoothing after T/2 steps. It is readily verified that the resulting smoothed random forcing error fields Σ1/2fwi 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.
  2. For each realization evaluate the tidal error field δu at the data locations xk and add measurement noise (scaled with the assumed data error) to compute a synthetic data vector di. 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.

  3. For each βi, use (21)–(23) to compute the inverse solutions ûi = u0 + δui (including elevations and transports).

  4. Using the actual synthetic solutions ui = u0 + δui 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 
    formula
    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 (M2, S2, K1, O1, N2, P1, K2, Q1) are modeled.

Fig. 2.

Topography of the model domain region. White area corresponds to water depth less than 100 m

Fig. 2.

Topography of the model domain region. White area corresponds to water depth less than 100 m

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.

Fig. 3.

Data locations used in study. TOPEX/Poseidon data sites shown as small open circles. Representer sites are shown with various filled black symbols, as described in text just prior to the callout of Fig. 9 

Fig. 3.

Data locations used in study. TOPEX/Poseidon data sites shown as small open circles. Representer sites are shown with various filled black symbols, as described in text just prior to the callout of Fig. 9 

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 M2 in the South China Sea is given in Fig. 4.

Fig. 4.

Example of an elevation representer for an observation of M2 for a data location in the China Sea indicated by the asterick. Amplitude contours are shown with solid lines, phase contours with dashed gray lines. The highest amplitude values (∼1.5 × 10−3) lie close to the representer location. The zero phase line goes through the data site since the representer must be real at this location

Fig. 4.

Example of an elevation representer for an observation of M2 for a data location in the China Sea indicated by the asterick. Amplitude contours are shown with solid lines, phase contours with dashed gray lines. The highest amplitude values (∼1.5 × 10−3) lie close to the representer location. The zero phase line goes through the data site since the representer must be real at this location

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.

Table 1.

Summary of computational resources needed for representer calculation for eight constituents in single precision

Summary of computational resources needed for representer calculation for eight constituents in single precision
Summary of computational resources needed for representer calculation for eight constituents in single precision

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 (3m + 4)Noc complex words, where Noc 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 M2 and K1 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 M2, and 40 cm for K1.

Fig. 5.

Elevation for the (a), (b) prior and inverse (c), (d) solution for the M2 (left) and K1 (right) tidal constituent. Amplitude (in cm) is shown with solid lines, and phases with dashed lines. The amplitude increment is 10 cm for M2 and 5 cm for K1. The phase increment is 90° for M2 and 45° for K1

Fig. 5.

Elevation for the (a), (b) prior and inverse (c), (d) solution for the M2 (left) and K1 (right) tidal constituent. Amplitude (in cm) is shown with solid lines, and phases with dashed lines. The amplitude increment is 10 cm for M2 and 5 cm for K1. The phase increment is 90° for M2 and 45° for K1

Figures 6 and 7 show in-phase and quadrature tidal currents for M2 and K1, 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 M2 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).

Fig. 6.

The M2 (a) in-phase and (b) quadrature currents. The inverse solution is shown in the larger images, with areas significantly different from the prior framed. Corresponding parts of the prior solution are shown in the lower rows of (a) and (b)

Fig. 6.

The M2 (a) in-phase and (b) quadrature currents. The inverse solution is shown in the larger images, with areas significantly different from the prior framed. Corresponding parts of the prior solution are shown in the lower rows of (a) and (b)

Major changes in the K1 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 M2 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.

Fig. 8.

The M2 along-track residuals for the real (a), (c) (left) and imaginary (b), (d) (right) inverse solution and TPXO.3 global solution. TOPEX/Poseidon tracks are shown with straight lines. Curves denote the values of the residuals—the bigger the residuals the bigger the distance between a track and a curve. Filled areas denote positive residuals and unfilled areas negative residuals

Fig. 8.

The M2 along-track residuals for the real (a), (c) (left) and imaginary (b), (d) (right) inverse solution and TPXO.3 global solution. TOPEX/Poseidon tracks are shown with straight lines. Curves denote the values of the residuals—the bigger the residuals the bigger the distance between a track and a curve. Filled areas denote positive residuals and unfilled areas negative residuals

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 M2 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).

Table 2.

Rms of along tracks residuals (cm) for different solutions and depth ranges

Rms of along tracks residuals (cm) for different solutions and depth ranges
Rms of along tracks residuals (cm) for different solutions and depth ranges

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.

Fig. 9.

Rms misfit vs number of representers used in the reduced basis solution. Curves for semidiurnal constituents are solid, and for diurnal constituents dotted: (a) shallow depth range, H < 500 m; (b) deeper depth range, 500 m < H < 2500 m

Fig. 9.

Rms misfit vs number of representers used in the reduced basis solution. Curves for semidiurnal constituents are solid, and for diurnal constituents dotted: (a) shallow depth range, H < 500 m; (b) deeper depth range, 500 m < H < 2500 m

Figures 10a,b show posterior errors for the M2 and K1 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.

Fig. 10.

Posterior error estimates. Dashed lines contour the area with error more than 1 cm. Outside of this area errors are less than 1 cm (over most of the open ocean). Solid contours are spaced with a 2-cm increment starting from 2 cm, and the area with error values exceeding 4 cm is filled with gray: (a) M2 and (b) K1

Fig. 10.

Posterior error estimates. Dashed lines contour the area with error more than 1 cm. Outside of this area errors are less than 1 cm (over most of the open ocean). Solid contours are spaced with a 2-cm increment starting from 2 cm, and the area with error values exceeding 4 cm is filled with gray: (a) M2 and (b) K1

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 105 grid nodes and 104 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.

Fig. 7.

The K1 (a) in-phase and (b) quadrature currents. The inverse solution is shown in the larger images, with areas significantly different from the prior framed. Corresponding parts of the prior solution are shown in the lower rows of (a) and (b)

Fig. 7.

The K1 (a) in-phase and (b) quadrature currents. The inverse solution is shown in the larger images, with areas significantly different from the prior framed. Corresponding parts of the prior solution are shown in the lower rows of (a) and (b)

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

REFERENCES
Anderson
,
E.
, and
and Coauthors
,
1992
:
Lapack User's Guide.
Society for Industrial and Applied Mathematics, 235 pp
.
Baines
,
P. G.
,
1982
:
On internal tide generation models.
Deep-Sea Res
,
29
,
307
338
.
Bell
,
T. H.
,
1975
:
Topographically induced internal waves in the open ocean.
J. Geophys. Res
,
80
,
320
327
.
Bennett
,
A. F.
,
1992
:
Inverse Methods in Physical Oceanography. Monographs on Mechanics and Applied Mathematics,.
Cambridge University Press, 346 pp
.
Bennett
,
A. F.
, and
P. C.
McIntosh
,
1982
:
Open ocean modeling as an inverse problem: Tidal theory.
J. Phys. Oceanogr
,
12
,
1004
1018
.
Cummins
,
P. F.
,
D.
Masson
, and
M. G. G.
Foreman
,
2000
:
Stratification and mean flow effects on diurnal tidal currents off Vancouver Island.
J. Phys. Oceanogr
,
30
,
15
30
.
Dushaw
,
B. D.
,
G. D.
Egbert
,
P. F.
Worcester
,
B. D.
Cornuelle
,
B. M.
Howe
, and
K.
Metzger
,
1997
:
A TOPEX/Poseidon global tidal model (TPXO.2) and barotropic tidal currents determined from long-range acoustic transmissions.
Progress in Oceanography, Vol. 40, Pergamon, 337–367
.
Egbert
,
G. D.
,
1997
:
Tidal data inversion: Interpolation and inference.
Progress in Oceanography, Vol. 40, Pergamon, 81–108
.
Egbert
,
G. D.
, and
A. F.
Bennett
,
1996
:
Data assimilation methods for ocean tides.
Modern Approaches to Data Assimilation in Ocean Modeling, P. Malonotte-Rizzoli, Ed., Elsevier Science, 147–179
.
Egbert
,
G. D.
, and
R. D.
Ray
,
2000
:
Significant dissipation of tidal energy in the deep ocean inferred from satellite altimeter data.
Nature
,
405
,
775
778
.
Egbert
,
G. D.
, and
R. D.
Ray
,
2001
:
Estimates of M2 tidal energy dissipation from TOPEX Posiedon altimeter data.
J. Geophys. Res., 106, 22 475–22 502
.
Egbert
,
G. D.
,
A. F.
Bennett
, and
M. G. G.
Foreman
,
1994
:
TOPEX/Poseidon tides estimated using a global inverse model.
J. Geophys. Res
,
99
,
24821
24852
.
Farrell
,
W. E.
,
1972
:
Deformation of the earth by surface loads.
Rev. Geophys. Space Phys
,
10
,
761
797
.
Golub
,
G. H.
, and
C. F.
Van Loan
,
1989
:
Matrix Computation.
2d ed. Johns Hopkins University Press, 642 pp
.
Hatayama
,
T.
,
T.
Awaji
, and
K.
Akimoto
,
1996
:
Tidal currents in the Indonesian Seas and their effect on transports and mixing.
J. Geophys. Res
,
101
,
12353
12373
.
Hendershott
,
M. C.
,
1972
:
The effects of solid earth deformation on global ocean tides.
Geophys. J. Roy. Astron. Soc
,
29
,
389
402
.
Hendershott
,
M. C.
,
1977
:
Numerical models of ocean tides.
Marine Modeling, E. Goldberg, Ed., The Sea, Vol. 6, John Wiley and Sons, 47–96
.
Kantha
,
L. H.
,
1995
:
Barotropic tides in the global oceans from a nonlinear tidal model assimilating altimetric tides. Part I: Model description and results.
J. Geophys. Res
,
100
,
25283
25308
.
Kivman
,
G. A.
,
1997
:
Assimilating data into open ocean tidal models.
Surv. Geophys
,
18
,
619
643
.
Koblinsky
,
C. J.
,
R. D.
Ray
,
B. D.
Beckeley
,
Y-M.
Wang
,
L.
Tsaoussi
,
A.
Brenner
, and
R.
Willimanson
,
1999
:
NASA ocean altimeter Pathfinder Project Report 1.
Data processing handbook. NASA/TM-1998-208605, 55 pp
.
Ledwell
,
J. L.
,
E. T.
Montgomery
,
K. L.
Polzin
,
L. C.
St. Laurent
,
R. W.
Schmitt
, and
J. M.
Toole
,
2000
:
Evidence for enhanced mixing over rough topography in the abyssal ocean.
Nature
,
403
,
179
182
.
Le Provost
,
C.
, and
A.
Poncet
,
1978
:
Finite element method for spectral modeling of tides.
Int. J. Numer. Methods Eng
,
12
,
853
871
.
Lyard
,
F. H.
,
1999
:
Data assimilation in a wave equation: A variational representer approach for the Grenoble tidal model.
J. Comput. Phys
,
149
,
1
31
.
Matsumoto
,
K.
,
M.
Ooe
,
T.
Sato
, and
J.
Segawa
,
1995
:
Ocean tide model obtained from TOPEX/Poseidon altimetry data.
J. Geophys. Res
,
100
,
25319
25330
.
Mazzega
,
P.
, and
M.
Berge
,
1994
:
Ocean tides in the Asian semienclosed seas from TOPEX/Poseidon.
J. Geophys. Res
,
99
,
24867
24881
.
McIntosh
,
P. C.
, and
A. F.
Bennett
,
1984
:
Open ocean modeling as an inverse problem: M2 tides in Bass Strait.
J. Phys. Oceanogr
,
14
,
601
614
.
Munk
,
W. H.
, and
C.
Wunsch
,
1998
:
Abyssal recipes II: Energetics of tidal and wind mixing.
Deep-Sea Res
,
45
,
1977
2010
.
National Geophysical Data Center
,
1992
:
Worldwide marine geophysical data.
GEODAS CD-ROM NOAA 92-MGG-02
.
Parker
,
R. L.
,
1994
:
Geophysical Inverse Theory.
Princeton University Press, 386 pp
.
Parker
,
R. L.
, and
L.
Shure
,
1982
:
Efficient modeling of the earth's magnetic field with harmonic splines.
Geophys. Res. Lett
,
9
,
812
815
.
Pingree
,
R. D.
,
1983
:
Spring tides and quadratic friction.
Deep-Sea Res
,
30
,
929
944
.
Polzin
,
K. L.
,
J. M.
Toole
,
J. R.
Ledwell
, and
R. W.
Schmitt
,
1997
:
Spatial variability of turbulent mixing in the abyssal ocean.
Science
,
276
,
93
96
.
Ray
,
R. D.
,
1998
:
Ocean self-attraction and loading in numerical tidal models.
Mar. Geod
,
21
,
181
192
.
Ray
,
R. D.
,
1999
:
A global ocean tide model from TOPEX/Poseidon altimetry.
NASA Rep. NASA/TM-1999-209478, 58 pp
.
Ray
,
R. D.
,
2001
:
Inversion of oceanic tidal currents from measured elevations.
J. Mar. Syst
,
28
,
1
18
.
Ray
,
R. D.
, and
B. V.
Sanchez
,
1989
:
Radial deformation of the earth by oceanic tidal loading.
NASA Tech. Memo. 1000743, 49 pp
.
Ray
,
R. D.
, and
G. T.
Mitchum
,
1997
:
Surface manifestation of internal tides in the deep ocean.
Progress in Oceanography, Vol. 40, Pergamon, 135–162
.
Ray
,
R. D.
,
R. J.
Eanes
,
G. D.
Egbert
, and
N. K.
Pavlis
,
2001
:
Error spectrum for the global M2 ocean tide.
Geophys. Res. Lett
,
28
,
21
25
.
Schwiderski
,
E. W.
,
1980
:
Ocean tides. Part II: A hydrodynamic interpolation model.
Mar. Geod
,
3
,
219
255
.
Shum
,
C. K.
, and
and Coauthors
,
1997
:
Accuracy assessment of recent ocean tide models.
J. Geophys. Res
,
102
,
25173
25194
.
Sjoberg
,
B.
, and
A.
Stigebrandt
,
1992
:
Computations of the geographical distribution of the energy flux to mixing processes via internal tides and the associated vertical circulation in the ocean.
Deep-Sea Res
,
39
,
269
291
.
Smith
,
W. H. F.
, and
D. T.
Sandwell
,
1997
:
Global seafloor topography from satellite altimetry and ship depth soundings.
Science
,
277
,
1956
1962
.
Yosida
,
K.
,
1980
:
Functional Analysis.
6 ed. Springer-Verlag, 500 pp
.

APPENDIX A

Numerical Implementation Details

We consider operators 𝗚 (gradient), 𝗗 (divergence), and 𝗖 (multiplication by Ω−1) as defined in section 3, Eq. (17). The nodes ζij, Uij, and Vij are defined as in Fig. 1, section 3. Depth Hij is defined on the ζ nodes. Let MU, MV, and Mζ be the water/land masks (0 for land nodes, 1 for ocean nodes). The node Uij or Vij 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

 
formula

Here a is the radius of the earth, g is gravitational acceleration, (θ, ϕ) are the latitude and longitude. The 2 × 2 matrices ΩUij and ΩVij are defined by (14) at the Uij and Vij nodes, and the following components of the inverse matrices are used in (A3):

 
formula

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:

 
𝗖: UijUij,  VijVij

and 𝗔 is identity on ζ boundary nodes:

 
𝗔: ζijζij.

It is easily verified from (18) that with these conventions, the solution on any boundary node (say Uij) will satisfy Uij = (fU)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

 
formula

the solution to 𝗔x = b satisfies

 
formula

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

 
formula

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 𝗔11y = b1 for y and 𝗔11zk = αk for zk, k = 1, m using the factorization of 𝗔11 (the factorization need be done only once). Collect zk, k = 1, m into the matrix 𝗭 = |z1, z2, … , zm|.

(b) Compute w = b2 − 𝗔21y and the square m × m matrix 𝗘 = 𝗔22 − 𝗔21𝗭.

(c) Solve the full m × m system 𝗘x2 = w.

(d) Compute u = b − 𝗔12x2.

(e) Solve 𝗔11x1 = 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

 
formula

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

 
formula

In the last step we have converted the real matrix 𝗕k back to complex via 𝗕k = 𝗕k1i𝗕k2, where 𝗕k = [𝗕k1 𝗕k2] is partitioned into two real submatrices, each with Nc columns. With the transformation of (C2) we have reduced the data vector to 2Nc components, discarding the part of the data that cannot be fit by any Nc constituent tidal model. Using the usual least squares (LS) theory, we estimate the data error variance from the mean-square residual magnitude σ̂2k = ‖dkk2/(Tk − 2Nc). Assuming that the original data errors ɛk are uncorrelated with variance σ2k, 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 σ̂2k to estimate the error variance for the 2Nc elements of k. For each data location xk this initial data reduction procedure yields a reduced data vector k, temporal sampling structure matrix 𝗕k, and data error variance σ̂2k. These can then be used along with the calculated representers to solve simultaneously for all Nc constituents as outlined in section 4. Note that if the time series at xj is shorter than 2Nc, 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 σ2k 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 dk, we can simply solve the triangular system k = 𝗕zk. The result will be the usual LS estimate of the real and imaginary parts of the Nc tidal harmonic constants. Error variances for the elements of zk are obtained in the usual way for the LS estimates from the diagonal elements of σ̂2k(𝗕Tk𝗕k)−1. These estimates of error variances for each constituent are then used to define the diagonal data error covariance Σ2e.

Next, we consider further details on matrix manipulations required for the multiple constituent inversion. Let Bk = |b1k b2k · · · bNck| denote the Nc columns of the complex matrices defined above. Define the large sparse blocked matrix

 
formula

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 Σ−1e𝗙, 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 NT eigenvalues. Typically NT can be approximately half of K with little effect on the final approximate inverse solution. With this reduction Σ−1e𝗙 is M × 2NcNT. Since M can be as large as 105, 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 Σ−1e𝗙 = 𝗪F𝗤F, where 𝗤F is an upper-triangular matrix of size 2NcNT × 2NcNT. Computation of 𝗤F and 𝗪TF 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 Σ−1e𝗙 = (𝗪F𝗪2)𝗦𝗤T = 𝗪𝗦𝗤T is the SVD of the large matrix Σ−1e𝗙. Note that the product 𝗪T = 𝗪T2𝗪TF [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