A near-global model for the sea surface expression of the baroclinic tide has been developed using exact-repeat mission altimetry. The methodology used differs in detail from other altimetry-based estimates of the open ocean baroclinic tide, but it leads to estimates that are broadly similar to previous results. It may be used for prediction of the baroclinic sea level anomaly at the frequencies of the main diurnal and semidiurnal tides , , , and , as well as the annual modulates of , denoted and . The tidal predictions are validated by computing variance reduction statistics using independent sea surface height data from the CryoSat-2 altimeter mission. Typical midocean baroclinic tidal signals range from a few millimeters to centimeters of elevation, corresponding to subsurface isopycnal displacements of tens of meters; however, in a few regions, larger signals are present, and it is found that the present model can explain more than 13-cm2 variance at some sites. The predicted tides are also validated by comparison with a database of hourly currents inferred from drogued surface drifters. The database is large enough to permit assessment of a simple model for scattering of the low-mode tide. Results indicate a scattering time scale of approximately 1 day, consistent with a priori estimates of time-variable refraction by the mesoscale circulation.
Satellite altimetry has enriched our understanding of ocean dynamics by providing a sustained and near-global view of mean sea level and mesoscale eddies during the last 25 years (Fu and Cazenave 2001). It is now widely used in ocean forecasting (Willis et al. 2010), and it is contributing to a broad range of research on ocean and climate processes (Łyszkowicz and Bernatowicz 2017). Studies of ocean tides have been invigorated by the unique datasets generated with altimetry, leading to improved tide predictions (Stammer et al. 2014) and quantitative maps of tidal dissipation (Egbert and Ray 2000, 2001).
The astronomical tidal forcing (i.e., the perturbation of near-Earth gravity caused by the sun and moon) does not vary significantly over the depth of the ocean (Doodson and Warburg 1941), so the work done on the ocean by the tidal forcing is input almost exclusively to the barotropic tide (Kelly 2016). Significant energy loss from the barotropic tide occurs in shallow water on continental shelves and in deep water at seafloor topography. In the latter case, the tidal flow disturbs isopycnal surfaces and creates baroclinic pressure gradients, which propagate as internal waves (Baines 1982). This process of barotropic tidal energy loss, leading to baroclinic tidal generation in the deep ocean, amounts to roughly a 1-TW rate of work (Egbert and Ray 2001, 2003; Egbert et al. 2004). Although the latter barotropic-to-baroclinic conversion accounts for only 1/3 of the tidal dissipation, the rate of work is similar to that done by wind on the ocean (Scott and Xu 2009), exclusive of the wind work involved in generating surface waves and stirring the mixed layer. The vertical potential energy flux driven by the dissipation of these energy inputs is connected to the thermohaline circulation, meridional heat flux, and other climatically relevant transport processes (Wunsch and Ferrari 2004). Once generated, the baroclinic tide can transport energy for thousands of kilometers, but the locations and mechanisms whereby it dissipates are not known adequately. Ocean circulation climate models are sensitive to the detailed spatial distribution of the dissipation (Melet et al. 2013), so there is ongoing concern with mapping it empirically.
Because tidal periodicities are aliased by satellite sampling, altimetry can only identify that component of the baroclinic tide that is phase locked with the astronomical tidal forcing. Hence, the baroclinic dissipation inferred from altimeter-derived maps is a combination of apparent dissipation due to the loss of coherence of the tide caused by time-variable modulations of the propagation medium (Park and Watts 2006; Rainville and Pinkel 2006) and irreversible dissipation due to processes such as wave–wave interactions (MacKinnon and Winters 2005; Ward and Dewar 2010; Wunsch 2017), shear-driven mixing (St. Laurent and Nash 2004), and wave breaking (Legg and Huijts 2006). These processes range from weakly to strongly nonlinear, and the degree to which they can be distinguished using ocean observations is not clear. Likewise, the extent to which turbulent transport resulting from these processes can be localized to the site of energy loss from the baroclinic tide is also not presently understood. The topic of baroclinic energetics is not directly addressed here; however, the discussion section indicates how a combination of sea surface height and tidal current observations might address these issues in the future.
The goal of the present work is to construct accurate maps of the baroclinic tide, useful for both tidal prediction and dynamical studies in the open ocean. Because in situ observations of baroclinic tidal currents are typically highly variable, it was a surprise when baroclinic tides were observed with altimetry (Ray and Mitchum 1996; Kantha and Tierney 1997; Ray and Cartwright 2001; Carrère et al. 2004). Although baroclinic tides are associated with subsurface isopycnal displacements up to 100 m (Alford et al. 2010), the baroclinic pressures are equivalent to just a few centimeters of ocean surface elevation, which makes them challenging to measure and reliably map. Nonetheless, some measure of success has been possible using data from long records of multimission satellite altimetry and data-fitting techniques that range from interpolation with radial basis functions (Ray and Zaron 2016) to plane wave fitting using theoretically predicted dispersion relations (Zhao et al. 2016) to Kalman filtering in the spatial domain (Dushaw 2015). The motivation to produce accurate baroclinic tidal predictions grows out of the desire to remove aliased tidal variability from a variety of ocean observations (Zaron and Ray 2017), which will be especially important for making use of observations from the anticipated Surface Water and Ocean Topography (SWOT) swath altimeter mission (Gaultier et al. 2016).
The rest of this paper is organized as follows. First, the altimeter data and harmonic analysis are briefly described, emphasizing minor innovations compared to previous approaches in the literature. Then, the spatial model and other details involved in estimating gridded tidal fields are described, along with a brief comparison with independent altimeter data to illustrate the use of the model for prediction of baroclinic tidal sea surface height (SSH). To infer energetics from the mapped SSH fields, some model for the tidal dynamics is needed, and this is explored by detailed comparison with a large database of surface currents observed with Lagrangian drifters. Finally, the results are discussed in the context of other studies, and the article concludes with a brief summary.
2. Altimetry data and harmonic analysis
The satellite altimeter data used in the present analysis are listed in Table 1, representing essentially all the exact-repeat altimeter mission data available during the 1992–2017 time period. By far the largest quantity of observations lie along the TOPEX/Jason reference orbit, but the other missions are essential for resolving the spatial structure of the baroclinic tide.
The path delay and geophysical corrections applied to the data are conventional and follow the Geophysical Data Record, version D (GDR-D), standard (Picot et al. 2012, 2014) with two minor innovations. The first innovation is that the barotropic ocean tide and Earth load tide are corrected using the Goddard/Grenoble Ocean Tide model, version 4.10c [GOT4.10c; an updated version of the model developed by Ray (1999)], which has been smoothly extrapolated by the author up to the coastline. The second innovation is that an estimate of the mesoscale sea level anomaly (SLA) is also subtracted from the SSH in order to remove as much nontidal signal as possible, prior to harmonic analysis. The rationale and impact of this correction have been discussed previously (Ray and Byrne 2010; Ray and Zaron 2016), so they will not be repeated here. As shown by Ray and Zaron (2011), the methodology used to estimate the mesoscale SLA from multisatellite altimetry (DUACS/Aviso Team 2014; Pujol et al. 2016) does not completely filter out the baroclinic tidal signals. Filtering the baroclinic tidal signals from the mesoscale SLA estimate prior to using it as a correction is necessary, and this procedure is detailed by Zaron and Ray (2018).
Aside from the minor changes just described, the data processing is essentially identical to that used previously for harmonic analysis of along-track altimetry. The altimeter data from missions with the same orbit ground tracks are assembled into time series at each point along the nominal ground track. Each time series is then harmonically analyzed using conventional methods (Cherniawsky et al. 2001; Carrère et al. 2004).
Some care is needed when choosing the tidal frequencies for mapping baroclinic tides. The first consideration is whether or not the aliased tidal frequencies can be determined with the given length of record for the given orbits. The second consideration is the extent of overlap or contamination by tidal frequencies that are not part of the analysis.
To address the first consideration, Table 2 lists the aliases and synodic periods for the tides that are mapped below: , , , , and the annual modulates of , denoted and (Huess and Andersen 2001). The synodic period is the time needed to accumulate a phase difference of 2π between signals at the given alias periods. For missions in the TOPEX/Jason reference orbit, the alias frequencies of these tides are unambiguously separable using time series longer than 3 years. For the Geosat Follow-On mission (G1A) in the Geosat orbit, the , , and frequencies can be accurately identified, but the and tidal aliases require a 12-yr time series to separate, and the and aliases essentially coincide. Thus, the estimates of , , and are inaccurate from G1A. The ERS/Envisat orbit is sun-synchronous, so cannot be determined by missions in this orbit. There is a near overlap of and , so these tides also cannot be determined accurately from these missions.
To address the second consideration, it is necessary to inspect the synodic periods of a larger set of tides. This has been done for the set , , , , , , , and , together with their annual modulates (not shown). For the missions in the TOPEX/Jason orbits, these additional frequency pairs extend the synodic period to about 6 years, which does not lead to problems with the present datasets. For the G1A orbit, the annual modulates of are found to nearly overlap with , which adds to noise in harmonic constants. Missions in the ERS/Envisat orbit are more problematic; separation of and its annual modulates from and its annual modulates requires almost 9 years of data (Andersen 1995), so the baroclinic tide is a source of noise here. Likewise, the – pair requires a nodal cycle to separate, so the presence of signals adds noise to . Overlaps of the tidal aliases and the annual and semiannual cycles also exist, but fortunately the annual and semiannual signals in the open ocean are captured by the mesoscale SLA maps and removed prior to harmonic analysis.
3. Spatial model and mapping
Previous efforts to map the baroclinic tides have used a variety of models to describe their spatial structure. For example, Ray and Zaron (2016) simply used an ad hoc radial basis function to smoothly interpolate harmonic constants between orbit ground tracks, yielding estimates of baroclinic tidal fields with a minimum of assumptions about their dynamics or spatial coherence. The consistency of the results with the predicted wavenumber dispersion relation for linear inertia–gravity waves supports other methodologies in which these dynamics are assumed. The maps of Dushaw (2015) directly use the dispersion relation for internal waves at the tidal frequencies, deriving the spatial coherence from the assumed dynamics. An even more constrained spatial model is the plane wave fitting used by Zhao et al. (2016), in which the spatial fields are assumed to comprise a small number of waves propagating in directions inferred from the data.
Experience with the plane wave fitting indicates that the baroclinic tidal fields closely obey linear dynamics (Ray and Cartwright 2001; Zhao et al. 2012; Zhao 2016); however, there is a tradeoff between bias and stability that depends on the complexity of the spatial model. For example, both empirical estimates and numerical models of baroclinic tides find a great deal of spatial structure and anisotropy, with waves organized into relatively narrow beams as the result of distributed sources and wave interference (Rainville et al. 2010). The unstructured signal model of Ray and Zaron (2016) is biased toward zero far from the data sites, and it is found that simply increasing the harmonic constants by 5%–20%, depending on location, improves the accuracy of the tidal predictions made with the model. Similarly, one would expect the highly structured plane wave model of Zhao et al. (2016) to be biased in a wave field composed of relatively narrow beams. The bias depends on the size of the fitting window, but one would expect it to be proportional to , where η is baroclinic tidal amplitude and y is a coordinate perpendicular to the local propagation direction.
Based on the above, a spatial model was hypothesized that represents the baroclinic wave field locally with a small number of propagating waves combined with a polynomial amplitude modulation. To make these ideas precise, let represent Cartesian coordinates on a locally defined tangent plane, and assume that the baroclinic tide can be represented with N spatially modulated plane waves, each with wavenumber modulus and direction , for . Assuming the amplitude envelope is modulated by a polynomial of order P, then the local spatial signal model for the complex amplitude of the baroclinic tide η is given by
where vector wavenumber is given by , and complex coefficients are found by weighted least squares fitting to the harmonically analyzed altimeter data. With this representation, the component of the wave field propagating parallel to is given by
and the antiparallel component is given by
Note that η and the model parameters are together a function of tidal frequency . When it is necessary to indicate this dependence, it will be shown using superscript notation [e.g., is the complex amplitude of the harmonic constant].
It remains to be stated how P, N, kn, and ϕn are determined. The procedure is explained here, but the reader may wish to note Fig. 1, which indicates the regions illustrated in detail in Figs. 2–4 referred to below. The order of the polynomial is simply chosen as P = 2, which eliminates the leading-order bias term of a plane wave fit. Determining N, the number of component waves, is done with a preliminary analysis in which and are determined through an incremental model-building exercise. The procedure is as follows. The empirical along-track harmonic constants are assembled in locally defined coordinates and , where L is the size of the two-dimensional data-fitting window. Within the fitting window, the harmonic constants are placed in square bins of size , averaging data from crossing tracks if necessary. The data within the grid are then regarded as the field of harmonic constants multiplied by the spatial sampling pattern of the altimeter ground tracks (Fig. 2a). The contents of this array are windowed (Fig. 2b), and the two-dimensional Fourier transform is taken. In essence, the resulting two-dimensional spectrum is the convolution of the baroclinic tide SSH with the antenna pattern of the ground tracks (Fig. 2c). In spite of the modulation by the antenna pattern, peaks in the spectrum are clearly identifiable (Fig. 2d). The two-dimensional spectrum is integrated azimuthally, and the peak wavenumber is used to assign (Fig. 2g). Then, the two-dimensional spectrum is integrated radially from to , and the azimuthal direction of the peak is assigned to (Fig. 2j). A simple plane wave fit is computed and subtracted from the data, and the process is repeated until an insignificant amount of variance is removed.
It is useful to examine examples of this procedure in different regions, and this is shown in Figs. 2–4 for three very different sites in the Pacific. Figure 2 is from a region northeast of the Hawaiian Ridge, where the wave field is dominated by a northbound mode-1 wave from the Ridge and a southbound wave from the Aleutians. The power spectrum of the windowed data (Fig. 2d) clearly shows the peaks associated with the northward- and southward-propagating beams, with the expected wavelength (170 km; Fig. 2g) and propagation directions (about ±70° from east; Fig. 2j). When the first wave is removed, the spectrum of the residual is dominated by the southward wave (Figs. 2e,h,k). The splitting of the northward peak after the first and second waves are removed (Figs. 2e,f) indicates that it is not well represented by a simple plane wave. Based on the shape of the spectral peak, it appears to be better represented by a radially spreading wave (Fig. 2d); however, this spatial model is not part of this preliminary exercise, which is only intended to identify . The units of the integrated spectra in the last two rows, mm2 (cpk)−1 and mm2 (rad)−1, where cpk is cycles per kilometer, allow the results to be compared with the data in the following two figures. Note that the two-dimensional spectra in the second row are log-scaled, showing three orders of magnitude, and the colors are normalized by the maximum value. Thus, as the residual gets smaller, the peaks stand out less above the noise floor (Figs. 2f,i,l).
The wave field in the subtropical western Pacific is primarily composed of three mode-1 waves (Figs. 3a,b,d), and the amplitudes of these waves (Figs. 3g–i) are larger than those described above near the Hawaiian Ridge. Note also that the noise floor of the radial wavenumber spectrum (e.g., Fig. 3i) is noticeably elevated, compared to the previous case (Fig. 2i). Presumably, this is related to the higher level of mesoscale kinetic energy in this region.
The wave field in the equatorial Pacific is much more directional than the previous examples (Figs. 4a,b,d). The wave fitting identifies two mode-1 waves and one mode-2 wave, all propagating to the south. Once again, notice how fitting with a single plane wave changes the directional distribution of variance (Fig. 4j vs Fig. 4k), suggesting that radial spreading is significant even within these small analysis windows.
The antenna patterns are generally similar in these examples since they primarily depend on the ground track spacing among the missions (Figs. 2c, 3c, 4c). Note that the Nyquist wavenumber, approximately 0.08 cpk, lies far outside the displayed range of wavenumbers. In fact, the antenna patterns are increasingly structured at large wavenumbers because of the high wavenumbers associated with the across-track sampling. Fortunately, the tidal fields contain so little variance at these small spatial scales that the leakage is not problematic for low wavenumbers, .
Table 3 lists the parameters for the spatial models for each of the tides considered. The analysis window of , , is smaller than that used for the other tides. A larger window, , is used for , , and because the along-track estimates of these tides are less accurate than , as discussed in section 2. The larger window, , is used for the and tides in order to resolve the longer wavelength of these diurnal tides, compared to the tide. Although the window is of size L, parameters in the model are determined by fitting the data with a weighting function, , so essentially just data from the middle third of the analysis window are used. As mentioned in the caption of Table 3, this function is also used to window the data prior to computing the two-dimensional power spectrum for determination of .
The procedure just described leads to a sequence of estimates for the dominant wavenumbers, modulus and direction , for , ordered according to the variance explained in the two-dimensional wavenumber domain. But how should the size of this expansion N be determined? Experimentation with the F test, in which the ratio of explained-to-prior variance is compared to that expected by chance (Jenkins and Watts 1968), and Aikake’s Information Criterion (Bozdogan 1987) found that both methods sometimes lead to spurious results, apparently due to occasional outliers. Instead, a simpler criterion was used. Namely, the expansion was truncated at when the wave removed less than 1.5-mm2 variance. The numeric value here was chosen to be approximately equal to the formal error estimate of the harmonic constants from the longest merged time series (TXA, J1A, J2A, J3A; Table 1). The value of N, which is typically in the range of two to five waves, is thus based on a subjective criterion designed to avoid overfitting the observations.
With the spatial model defined as above, the mapping proceeds by dividing the ocean into patches of size in a local tangent plane centered on latitude and longitude coordinates . The patches lie on a regular overlapping grid of latitudes, , such that , where is the mean radius of Earth, and determines the extent of overlap. The longitude grid is also equidistant between the local tangent planes , where ; note that and depend on l, but this dependence is suppressed in the notation for readability.
Previous published maps of the baroclinic tides have utilized along-track high-pass filtering of the data in order to suppress errors at wavelengths longer than 500 km, but this leads to a nonisotropic antenna response and filtering of east–west-propagating waves (Ray and Zaron 2016). To overcome this problem, the present approach estimates the model parameters by fitting the along-track sea surface slope, rather than SSH. This reduces the influence of long-wavelength errors in the data, but because the same operation is applied to both the input data and the signal model, there is no loss of sensitivity to wavenumbers oriented in the east–west directions.
Finally, the complex coefficients in Eq. (1) are determined by conventional weighted least squares within each analysis window. The weights used are the inverse of the squared standard error estimates from the along-track harmonic analysis (Cherniawsky et al. 2001).
Once the model coefficients are found for each patch, the estimated tidal fields are gridded on a regular latitude–longitude grid at a resolution of by weighted averaging of the overlapping patches. The averaging kernel is a radial basis function , so the resulting field is essentially continuously differentiable at the edges of each patch.
One final step is involved in preparing a high-resolution grid suitable for making tidal predictions, which is masking off regions where the estimate is thought to be inaccurate. This is done using the formal error estimate of the harmonic constants determined from altimetry, averaged over 500 km. The mask is set to zero where the mean standard error is greater than , a value that was subjectively determined. Additional criteria that result in a grid cell being masked off are the following: 1) fewer than 1250 data points used in the determination, 2) water depth less than 500 m, based on the General Bathymetric Chart of the Oceans (GEBCO) bathymetry (Weatherall et al. 2015), 3) less than 12-km distance to land, and 4) poleward of 60° latitude. As a final step, the discrete-valued mask is convolved with a compactly supported, twice-continuously differentiable function (Wendland 1995) in order to smooth the mask over 3°. In general, the mask delimits regions where the mapped field appears to be spurious; however, it would be advantageous to optimize the mask using more objective criteria in the future.
4. Assessment of baroclinic tide estimates
Figure 5 shows the , , , and tides obtained with the approach described above, plotting the elevation in phase with the Greenwich phase of the astronomical tidal potential (Simon 2013; detailed plots of the tides are available in the online supplemental material). The estimates obtained for the , , , and baroclinic tides are superficially similar to those shown in previous works (Dushaw 2015; Zhao et al. 2016; Ray and Zaron 2016; Zhao 2017), but they differ in quantitative detail. The most conspicuous difference is the better representation of the tide in the western Pacific, where the spatial aliasing of the tidal wavelength on the satellite ground tracks caused it to be reduced in previous efforts that utilized along-track spatial smoothing. In addition, the signal model of the present approach appears to admit more small-scale detail than previous estimates, although a dedicated intercomparison effort is still ongoing (L. Carrère 2018, personal communication; L. Carrère et al. 2018, meeting presentation). The detail visible in the field (Fig. 5b) is considerably reduced, compared to that of . This occurs because the map uses less data for , compared with , but also the amplitude of is lower and closer to the noise level. The way the spatial model is constructed essentially has a small-signal cutoff to avoid overfitting the data.
Maps of the annual modulates of in Fig. 6 are a new component of this work. The baroclinic and tides are generally too small to estimate reliably, so a larger fitting window has been used, and the result is heavily weighted toward their values from missions in the TOPEX/Poseidon reference orbit. The Arabian Sea, the region between Seychelles and Madagascar, and the region offshore of the Amazon River plume are locations where seasonal modulations are detectable (Fig. 6a). Seasonal modulation of the internal tides in the South China Sea have been studied previously and attributed to the seasonal cycle of upper-ocean stratification (Fig. 6b; Jan et al. 2008). The present maps essentially provide a regional view of these changes, which are present throughout the western Pacific. Note that baroclinic and barotropic tidal seasonality has been identified previously using altimeter data (Müller et al. 2014).
A comparison of the present estimates of with a similar estimate published by Ray and Zaron (2016), also at (1/20)° resolution, is shown in Fig. 7. The present estimate is denoted “high-resolution empirical tide” (HRET), and the Ray and Zaron (2016) estimate is denoted “internal tide” (IT). The difference of the in-phase components, HRET minus IT, is smaller than a few millimeters over much of the ocean, but differences exceeding a centimeter occur in the western Pacific and in a few other regions where amplitudes are large and wavenumbers are zonal. The differences display a distinct pattern of satellite ground tracks and indicate that the estimates are in best agreement along the tracks (not shown); however, in regions where the wavenumber is zonal, the difference is not a random error, but it consists of propagating waves. This difference for zonal waves is thought to be caused by the along-track filtering used by Ray and Zaron (2016) to reduce the influence of long-wavelength errors, but it also tends to filter the waves oriented perpendicular to the satellite ground tracks. The present approach is based on fitting a model of sea surface slope, so no explicit along-track filtering is necessary.
Figure 8 illustrates the usefulness of the separate tide models for dealiasing tides in independent altimeter data. The variance reduction of the CryoSat-2 altimeter SSH measurements is plotted using data from 2012 to 2018, averaging within 2° latitude–longitude bins (more than 10 000 measurements per bin). Positive values (red) indicate the degree to which the predicted tides remove aliased tidal variability, while negative values (blue) indicate that noise is added by the tide model. The maximum variance explained within any bin is 13 cm2 for , but the mean variance is only 0.25 cm2. There are a few locations, particularly to the south of the Kuroshio, where the explained variance is negative, but this is a region where fewer CryoSat-2 data are available due to changes in its data collection mask, so the significance of these variance estimates is reduced. explains a much smaller amount of variance, essentially all within the ±30° latitude range displayed. The model for explains a maximum of 4-cm2 variance, almost all in the seas of the western Pacific.
The model is less successful in explaining variance at the , , and frequencies (Figs. 8d–f). These are smaller tides, but there are more regions where the tidal correction fails to reduce the variance. Nonetheless, the fields indicate a few regions where these tides are large enough that they might be considered for use as tidal corrections, depending on the specific application. The variance reduction from the total of the tidal corrections is dominated by the and components (Fig. 8g).
5. Baroclinic tidal dynamics
The previous section focused on the sea surface height expression of baroclinic tides. Potentially more insight into the dynamics can be obtained by studying the baroclinic tidal currents. Let represent the horizontal current vector at the ocean surface. The instantaneous tidal currents ought to be related to the surface elevation through the equations
where is a vector that is a nonlinear function of both tidal currents and nontidal currents, with the latter denoted by . The relationship between the mean, phase-locked tidal currents and the surface elevation is, in principle, more complicated because it involves the projection of the above dynamics onto particular tidal frequencies:
where it is understood that and now refer to complex-valued fields associated with the jth tidal frequency , and is analogous to the divergence of a Reynolds stress.
The physical effects represented by can be thought of as tidal self-interactions, such as shear-driven mixing (St. Laurent and Nash 2004), scattering by topography (Johnston et al. 2003), or sub- and superharmonic generation (MacKinnon and Winters 2005; Wunsch 2017); and tidal–mean flow interactions, such as time-dependent refraction (Rainville and Pinkel 2006; Park and Watts 2006) or directional scattering by geostrophic modes (Ward and Dewar 2010). For the essentially low-mode description of the tidal elevation, which can be inferred from altimetry, it is hypothesized that time dependence of the propagation medium is the dominant physical process, and it can be approximated by the linear relationship
where is a damping time scale. One can estimate from the effective diffusivity of the nontidal flow , where and are the root-mean-square phase speed perturbation and its decorrelation scale, respectively. Then, , where is the wavenumber of internal tide. Plausible estimates for range from 0.05 to 0.1 m s−1 (Zaron and Egbert 2014; Buijsman et al. 2016), with a correlation scale of 100–400 km (Zaron and Egbert 2014). Assuming a midlatitude value of for a mode-1 baroclinic semidiurnal tide, the value of ranges from 8 × 10−6 to 6 × 10−5 s−1, a range of 5%–50% of the frequency. For , one would expect to be about a factor of 4 smaller because of the approximately doubled wavelength.
A test of the hypothesized dynamics [Eqs. (5) and (6)] has been conducted by predicting baroclinic tidal currents and comparing them with currents inferred from surface drifters. The dataset consists of 96 million hourly current vectors from 12 000 drogued drifters, collected from 1995 to 2015 as a part of the NOAA Global Drifter Program (Elipot et al. 2016). Observed currents are compared with predicted tidal currents, and the variance reduction is used as a measure of the goodness of fit. Figure 9 shows the explained variance as a map, averaged within 2.5° bins, when no damping is assumed . For the largest and most accurately determined tides and , the model explains a positive amount of vector current variance almost everywhere.
Several points of interest can be noted from Fig. 9. Comparison of the observed root-mean-square surface speed (Fig. 9a) with either of the root-mean-square predicted speeds (Figs. 9b,d) indicates that the predicted tidal currents are generally a small fraction of the observed currents; however, in a few areas, such as in Luzon Strait, near Seychelles, near New Zealand, and off the North American west coast, the tidal currents may comprise 20% or more of the total. Also, since the tidal currents are related to the gradient of sea surface height, the currents look quite different from the surface elevation (Fig. 9b vs Fig. 5a). Close inspection of the explained current variance does highlight a few sites where the tide model has problems. The model adds variance (negative explained variance) at a few spots in the North Pacific and at several other locations near the coast (Fig. 9c; regions colored blue). The Gulf of Mexico is a region where the predictions are not accurate (Fig. 9e).
The GDP dataset is large enough that it can discriminate between small adjustments to the dynamics. For example, the latitude-dependent acceleration of gravity g, which varies by about 0.5% from pole to equator, is used in Eq. (5) (Moritz 2000). If a constant nominal value is used instead, , the area-averaged explained variance is reduced slightly.
Figure 10 shows the explained variance as a function of for the and tides. The predictions for the other tides are not accurate enough to usefully constrain the damping time scale. The explained variance is maximized for and (Fig. 10a), values that are within the range proposed above. As is evident from Fig. 9, the geographic distribution of and currents is very different, so these variance-maximizing values of are measuring different physical locations. If, instead, the explained variance is restricted to the latitude range that includes Luzon Strait, where both and are relatively large, the optimal values are and (Fig. 10b), approximately , as predicted. The ratio of these coefficients varies somewhat when averaging over different regions when other conditional averages are used; however, the general property of has been observed in every case examined.
The maps of the low-mode , , and baroclinic tides presented here appear to be an incremental refinement of other published estimates (Zhao et al. 2016; Ray and Zaron 2016; Zhao 2017). Compared to these previously published models, the HRET model involves small changes in the signal model, fitting the model to sea surface slope rather than height, a slight increase in the quantity of data, and improved estimation and removal of nontidal variability prior to harmonic analysis. The mapping methodology was inspired by the plane wave fitting approach implemented by Zhao et al. (2016), which it sought to generalize and improve.
It is interesting to compare the present approach with one that uses the dynamics directly (e.g., Zaron et al. 2009). For the tide, the present approach estimates a maximum of 144 parameters per (250 km)2 analysis window (6 wavenumbers × 2 modulus and direction × 6 spatial polynomial coefficients × 2 in-phase and quadrature parts). Typically, this involves about 3600 harmonic constants, or about 25 data points per parameter to be estimated. If only three wavenumbers were identified in the model-building stage, then there would be 50 data points per parameter. This differs from the dynamics-based approach, where the number of parameters to be estimated is determined by the model’s spatial resolution and the assumed decorrelation scale of model errors. If the in-phase and quadrature components of the errors in the horizontal momentum equations were to be estimated, and a correlation scale of 50 km were assumed, there would be at least 200 parameters to estimate. This is roughly twice as many parameters as with the present approach and easily 3–6 times as many in regions with relatively few waves. Thus, the highly structured signal model using relatively few degrees of freedom, with the vector wavenumber set by the preliminary model building (a nonlinear estimator), seems to be an advantage, compared to more general dynamics-based approaches (L. Carrère et al. 2018, meeting presentation).
Nonetheless, it is clear that the present approach has limitations that will cause it to lose utility near abrupt topography or where the baroclinic wave field deviates from the rudimentary signal model. Figure 11 illustrates the harmonic constants in a region of the western Pacific where the Philippine Sea meets the East China Sea. The boundary between the baroclinic waves (in the deep water) and their absence (on the continental shelf) is apparent. Capturing this spatial structure with the signal model of Eq. (1) is not possible, at least within analysis windows containing a sufficient quantity of data. Instead, it seems likely that subsequent improvements will result from using a more dynamically constrained approach, but with a highly structured error covariance model to reduce the number of parameters involved.
The damping time scale of 105 s (or 1.2 days) estimated in section 5 agrees with the residence time of 1–1.5 days estimated by Zhao et al. (2016), based on their mode-1 tide map and the rate of energy input by the barotropic tide (Egbert and Ray 2001, 2003). Of course, the present estimate is based on completely independent data and methodology, so it provides a check on the previous energy budget summaries (Wunsch and Ferrari 2004; Garrett and Kunze 2007).
The values of estimated above were obtained by maximizing a goodness-of-fit metric, which implicitly emphasizes those spatial regions with the largest baroclinic tidal kinetic energy. Thus, the area average damping time scale of the coherent tide could be greater or less than the inferred value, and it is difficult to place confidence limits on it, or map its spatial structure, without a detailed consideration of the physical mechanisms it represents. Such an analysis will be the subject of future studies.
The form of the damping, , was justified as a model for the loss of energy from the phase-locked tide, but it may alternately be regarded as an energy source for the non-phase-locked tide. The non-phase-locked tide obeys an energy equation, which could be written as
where is the group velocity, and and are the non-phase-locked and phase-locked wave energies, respectively, associated with the jth tidal frequency. This expression can be integrated over the deep ocean, bounded by the continental margins, to obtain
where is the reflection coefficient for the low-mode baroclinic tide at the continental margin, about 0.4 (Kelly et al. 2013); P is the perimeter of the deep ocean; and A is the surface area of the ocean. At this level of approximation, the quantities and are regarded as area averages, and is a damping time for the non-phase-locked time.
One expects the ratio to be some multiple of the reciprocal of the radius of Earth , and is a measure of the propagation distance of the phase-locked tide, only about using a mean group speed of and . Even with a generous estimate for the perimeter P, across which the tides are reflected or dissipate, the ratio does not seem much larger than about , so the first term in the above expression is apparently smaller than 1/8. If the ratio is to be smaller than 1, then the non-phase-locked tide must be rapidly damped, with the same size or larger than .
It is hard to reconcile the estimate with known mechanisms of dissipation for the low-mode tide. It is possible that estimates of from altimetry are biased low (Zaron 2015, 2017). They are obtained from analysis of variance, which means they are implicitly weighted toward areas with the largest signals, presumably the generation sites where the non-phase-locked tide would be smaller than average. Estimates from moorings are not inconsistent with a larger value, (Alford and Zhao 2007), but even this is not large enough to constrain . Perhaps the regime described by Weisberg et al. (1987) is more typical.
Alternately, it is possible that the same factors causing the loss of coherence of the internal tide also result in the formation of caustics where the higher amplitudes and nonlinearity could lead to a rapid transfer of energy into the broadband internal wave spectrum (Zhao and D’Asaro 2011; Dunphy and Lamb 2014). Laboratory studies with coherent internal wave sources find that the rapid transfer of energy out of the phase-locked waves is enabled by lateral inhomogeneity in the wave field (Bordes et al. 2012).
A new series of models for the phase-locked component of the low-mode baroclinic , , , and tides and the annual modulations of the tide have been developed. The models differ from previous efforts in minor respects. While development and intercomparison is an ongoing exercise, it appears that the present results are slightly more accurate than other published and unpublished models (L. Carrère et al. 2018, meeting presentation). For example, the root-mean-square variance reduction of the present solution exceeds that of Ray and Zaron (2016) by only about 2 mm globally; however, larger differences exceeding 2.5 cm are present at specific sites (L. Carrère 2018, personal communication). Further improvements in satellite altimetry and processing techniques, and innovations in mapping techniques, will certainly lead to further increases in accuracy in the future.
The primary purpose of this manuscript is to document the mapping technique and validate the tide models using a large surface current drifter dataset. The latter leads to a new estimate for the scattering rate of the phase-locked tide, with implications for the generation and dissipation of the non-phase-locked tide.
The present model should be useful for removing aliased tidal signals from satellite altimeter measurements and in situ measurements of various kinds (Zaron and Ray 2017). The former should facilitate more accurate mapping of mesoscale sea level anomalies (Fu et al. 2010) and the identification of barotropic tidal energetics within small regions (Zaron and Egbert 2006). Tandem studies of barotropic and baroclinic tidal energetics, based on the present models, can be expected to lead to more rigorous bounds on the energetics than presented in section 6 and should be considered in the future.
The along-track altimetry data were obtained using the Radar Altimeter Database System (RADS; http://rads.tudelft.nl/rads/rads.shtml), which is an effort of the Department of Earth Observation and Space Systems at TU Delft and the NOAA Laboratory for Satellite Altimetry. The surface drifter data were collected and distributed via the NOAA Global Drifter Program (http://www.aoml.noaa.gov/phod/dac/gdp_information.php). Computations were performed using the Portland Institute for Computational Science cluster computer, Coeus, acquired with support from NSF Award DMS1624776 and ARO Award W911NF-16-1-0307. Further support was provided by the National Geospatial-Intelligence Agency Award HM0177-13-1-0008 and NASA Awards NNX16AH88G and NNX17AJ35G.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-18-0127.s1.