Interactive comment on “ A harmonic projection and least – squares method for quantifying Kelvin wave activity

while Referee #1 asks that your method should be “at least applied to the equatorial Pacific and compared with estimates from Boulanger and Menkes (1995, 1999)”, Referee #2 takes this point further, disagreeing with the assertion that you are building on that work, while offering considerable detail in challenging a number of the assumptions made in constructing your method. Referee #1 has also read, and does support, these comments.


Introduction
The quantification of ocean variability associated with equatorial long waves is a topic of great importance for understanding the tropical ocean and its role in climate.Since the advent of satellite 25 altimetry, the surface manifestations of these waves and the wind forcing driving them have been tracked in datasets that now comprise over 20 years of continuous global coverage (e.g., Delcroix et al., 1994;Susanto et al., 1998;Boulanger and Menkes, 1999;Drushka et al., 2010).However, to use these observations to better understand the behavior of these planetary waves and their relationship to climate variability, analysis techniques are needed that target the specific signatures of Kelvin 30 and Rossby waves in satellite observations.In particular, the present study was motivated by a need to quantify the relative presence of upwelling vs. downwelling Kelvin waves in the equatorial Indian Ocean and along the coasts of Sumatra and Java, where they are influential in the evolution of Indian Ocean Dipole events (Delman et al., submitted).
A variety of techniques have been employed to quantify equatorial long wave activity from satel-35 lite observations; these range from the application of sophisticated data assimilation techniques to meridional projections of sea level anomaly (SLA) data.The data assimilation approaches generally use a linear wave-propagation model, along with Kalman filters (e.g., Miller and Cane, 1989;Fu et al., 1993) or adjoints (e.g., Thacker and Long, 1988;Long and Thacker, 1989a, b) to incorporate observations.These techniques are particularly useful for cases where observations are sparse 40 and error-prone, as is often the case for in-situ measurements, and also during the earlier years of satellite observations when spatial resolution was low (e.g., Geosat).As the spatial and temporal coverage of altimeter-derived remote sensing data increased, it was conceivable to estimate Kelvin and Rossby wave activity using solely meridional projections of SLA data, or a combination of SLA and current observations.Cane and Sarachik (1981) showed that vectors containing SLA and sur-45 face current profiles associated with a given vertical Kelvin wave mode and its associated meridional Rossby wave modes are orthogonal; this orthogonality provided the basis for an equatorial wave decomposition in numerous studies (e.g., Delcroix et al., 1994;Yuan et al., 2004;Yuan and Liu, 2009).Boulanger andMenkes (1995, 1999), BM9599 hereafter, also carried out a decomposition using only meridional projections of SLA data that were reasonably consistent with projections derived 50 from in-situ moorings.However, the decomposition of Kelvin and Rossby wave modes based on meridional projections of SLA alone are not orthogonal, and as Yuan et al. (2004) notes, this necessitates the inversion of an ill-conditioned matrix.An alternative approach used complex EOFs of SLA to separate Rossby and Kelvin wave signals in the equatorial Pacific (Susanto et al., 1998); one limitation of this method is that complex EOFs by definition constrain the along-waveguide 55 and across-waveguide length scales of the waves, while shallow-water theory only constrains the across-waveguide length scale.
Here we build on the methodology of BM9599, by using the approximate phase speeds as well as cross-waveguide profiles to isolate the Kelvin wave signal.Starting with the SLA meridional time, followed by a least-squares fit to optimize the non-orthogonal projection coefficients.The result is a Kelvin wave coefficient K that approximates Kelvin wave generation and dissipation along the waveguide, and can be used to track coastal as well as equatorial Kelvin waves.The method as described is focused only on an accurate representation of Kelvin (not Rossby) wave activity, though an extension of these techniques might enable a comprehensive decomposition of equatorial waves 65 (as discussed in Section 4).The paper is structured as follows: Section 2 describes the satellite data used, and the harmonic projection and least-squares method that results in the Kelvin wave coefficient K. Section 3 estimates errors associated with the computation of K using a Monte Carlo simulation, and discusses qualitative and quantitative analyses of satellite observations to assess how effectively K describes Kelvin wave activity along the Indian Ocean waveguide.Section 4 70 summarizes the strengths and weaknesses of the method, and considers the possibility of extending the method to quantify Rossby wave activity.

Data
Our methodology quantifies Kelvin wave activity using AVISO Ssalto/Duacs alongtrack SLA data, 75 specifically those from the TOPEX/Poseidon, Jason-1, and Jason-2 satellites.These satellites repeat their orbit over a given track approximately every 10 days, and the data have near-continuous coverage from September 1992 to December 2013.The reason for using alongtrack as opposed to the frequently-used gridded product is the increased spatial resolution in the along-track direction, ∼1/10 • compared to 1/3 • for gridded data.One of the advantages of this method is its utility 80 for tracking waves in their transition from equatorial to coastal Kelvin waves.However, quantifying coastal Kelvin waves requires higher spatial resolution, as the baroclinic radius of deformation shrinks from ∼400 km at the equator to ∼100 km at 10 • S. The disadvantage of using the alongtrack data is the large spacing between tracks in the zonal/alongshore direction (∼ 300 km along the coast), but the spacing is still small relative to the along-waveguide length scale of Kelvin waves 85 near the equator, typically >1000 km.
Due to the anisotropy of equatorial-coastal long waves, the offset angle between satellite tracks and meridional cross-sections at the equator is likewise considered to be negligible, and both ascending and descending tracks are used in the analysis.Along the Sumatra and Java coasts, only ascending (SW-NE oriented) tracks are used in Kelvin/Rossby wave projections to best approxi-90 mate a cross-shore profile.For computational expediency in the least-squares part of the solution, the method was applied to overlapping two-year subsets of the full data record, with each subset overlapping with the next one by a year.The results from each subset were then patched together using a tapered weighted averaging in the overlapping year to create a continuous field of K values Ocean Sci.Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -1, 2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.
for the 21-year period of record (i.e., with 20 subsets patched together).For comparison purposes 95 and to present clear visual snapshots of variability in the Indian Ocean basin, gridded maps of SLA (MSLA) (Ducet et al., 2000) were also used to generate some of the figures in this paper.

Kelvin wave y-projections
The first step in the computation of the Kelvin wave coefficient K is to calculate the projection of the SLA data onto a meridional or cross-shore profile of a baroclinic Kelvin wave based on 100 linear shallow-water theory (e.g., Gill and Clarke, 1974;McCreary, 1981).We refer to this as the y-projection; for an equatorial Kelvin wave it is the same Gaussian profile given in Appendix A2 of Boulanger and Menkes (1995), but our analysis also considers coastal Kelvin waves for which the wave profile transitions to a decaying exponential away from the equator.For an equatorial-coastal Kelvin wave the profile is where y is the perpendicular distance relative to the equator or the coastline, h 0 is the amplitude (i.e., peak value) of the wave, f 0 is the Coriolis parameter at the latitude where the profile intersects the coast, and φ is the angle of orientation of the coast relative to the east-west axis (f 0 = 0 and φ = 0 for equatorial Kelvin waves).The sign in front of (f 0 /c)y for coastal Kelvin waves is chosen such that the term is always negative.As our focus here is on Indian Ocean Kelvin waves that are 110 deflected to the south of the equator, y is negative and decreasing away from the coast, and f 0 < 0, so the sign is negative.The value of c for the meridional/cross-shore profiles in this analysis was taken to be 2.5 m s −1 .This value of c lies between the first-and second-mode baroclinic phase speeds for Kelvin waves in the region, as these two modes account for most Kelvin waves observed in Indian Ocean SLA (Drushka et al., 2010).However, using c = 2.0 m s −1 or 3.0 m s −1 does not produce a 115 substantially different result.
Applied to the altimetry data, the Kelvin wave y-projection is given by for equatorial Kelvin waves and for coastal Kelvin waves south of the equator, where h SLA is the alongtrack altimetry profile, and r is the radius for the profile projection.The overbar indicates the meridional a = 1/(2r) cross-shore a = 1/(r) 0 −r a dy mean (for equatorial and coastal waves respectively) of the profile a over the range being integrated.For the r value, we used 5 • of latitude for equatorial Kelvin waves; r was then tapered to a distance equivalent to 3 • of latitude along the coasts of Java and Nusa Tenggara to account for the smaller radius of deformation.We note that K y is an integrated measure of the sea level displacement; this type of measure is a more consistent indicator of Kelvin wave activity in 125 the equatorial-coastal transition than peak amplitude, since without substantial dissipation, the peak amplitude of the wave tends to increase poleward as the radius of deformation decreases (Figure 1).

Projection using harmonic basis functions in x and t
After the Kelvin wave y-projections K y are computed, the next step in our approach is to project K y onto two-dimensional Fourier basis functions in along-waveguide distance x and time t.One 130 method of separating these components is to assume that a vector b consisting of the alongtrack Kelvin wave projections K y can be explained as a linear combination of two-dimensional Fourier basis functions where the columns of A are the basis functions A cos m,n = cos [2π (k m x − f n t)] and A sin m,n = sin [2π (k m x − f n t)] and x and t are along-waveguide distance and time respectively; the Fourier coefficients to be solved 135 for are contained in the vector m.
Basis functions A m,n that propagate from one side of the basin to the other at constant amplitude are most effective at representing Kelvin waves that similarly propagate across the basin with little change in amplitude.Kelvin waves that are forced and dissipate within the domain, especially with the low wavenumbers common to Kelvin waves, may have some of their energy aliased into 140 westward-propagating signals.To resolve this aliasing issue, we introduce an additional tapering parameter s to the basis functions (Figure 2).The basis functions A m,n,s take the form The tapering location x s is varied at intervals of ∆x = 600 km throughout the span of the waveguide, corresponding to the shortest wavelengths resolved along the coastal part of the waveguide (along the equator the effective Nyquist wavenumber is higher with more satellite tracks used).For s = 1, where D is the identity operator for s = 1 projections, and the finite difference operator for s > 1 projections.The column vector w can be used to weight the elements of Dm nP relative to one another.In this case setting w to all ones was found to be sufficient, though accuracy may be gained in certain areas by adjusting the weighting vector; we speculate about one such case in Section 4.
With w n = 1 for all n, we miminize the misfit of Dm nP and the size of m T m using the cost function Setting λ = 0.1 was found heuristically to produce the most credible reconstructions of the K y values in b, while reducing noise at the highest wavenumbers.The coefficient vector m is then given by Finally, all coefficients of m that correspond to westward-propagating basis functions, i. and K values for the same period (Figure 3).The K y and K values are calculated from the alongtrack SLA data at points where the satellite tracks cross the waveguide, hence these values are presented as points in Fig. 3b-c, whereas the raw SLA data from the gridded product are contoured in Fig. 3a.During the May-July period, the predominant feature in the raw SLA (Fig. 3a) is an eastward- propagating patch of elevated positive SLA, indicative of a downwelling Kelvin wave.However, the Kelvin wave y-projection (Fig. 3b) shows that this downwelling wave was both preceded and followed by upwelling Kelvin waves, both of which are much more evident in the Kelvin wave projection (Fig. 3b) than in the raw SLA (Fig. 3a).The y-projection still contains a number of and substituting into (1).A comparison of the reconstructed h K with gridded maps of SLA over a two-month period in 1997 (Figure 4) confirms that the Kelvin wave reconstruction is broadly consistent with the Kelvin wave activity suggested by the gridded SLA field, but also highlights some key differences.In late May and early June, an elevated SLA field persists in the eastern equato-210 rial Indian Ocean (Fig. 4b-c), while the reconstructed h K indicates that the Kelvin wave activity is changing sign from positive to negative there (Fig. 4g-h).In late June and early July, reconstructed h K indicates that upwelling Kelvin wave activity is being generated from approximately 60 • E eastward (Fig. 4i-j), while SLA is only substantially depressed east of 90 • E (Fig. 4d-e).The discrepancy during this latter period is likely accounted for by the influence of downwelling Rossby waves on 215 the SLA field in the central Indian Ocean; these waves have positive SLA maxima near the north and south radii of deformation, and would also elevate SLA (to a lesser extent) at the equator.Therefore, this implies that in the raw SLA field in early July 1997 (Fig. 4e) the upwelling Kelvin wave still present in the central equatorial Indian Ocean would not be apparent; this has potential implications for understanding the timing of the upwelling wave and where it was forced.220

Monte Carlo-based error estimates
In order to place uncertainty bounds on the method's capacity to remove westward-propagating wave activity from the Kelvin wave estimate, we carried out a Monte Carlo simulation.In this way the method could be applied to propagating waves whose amplitudes and K values were known a priori.
We created 100 years of randomly-generated basis function coefficients m, using Cholesky factor- Once the artificial wave field was constructed, the harmonic projection and least-squares method described in sections 2.3-2.4 was applied to the artificial K y field, and the K values derived from the basis function coefficients known a priori and deduced from the method were compared.An example of this for a given simulated year is shown in Figure 5; the artificially-generated K y field contains 240 signals propagating in both directions, though for most of the year the westward-propagating Rossby waves appear to predominate (Fig. 5a).However, a consideration of the Kelvin wave signal K a priori in isolation (Fig. 5b) reveals that in addition to the very strong downwelling wave early in the year, a series of weak and moderate Kelvin waves propagate throughout the year.Many of these weaker waves are unidentifiable in the K y field with the Rossby wave signals superimposed (Fig. 5a).However, the 245 reconstructed Kelvin wave signal K reconst , computed by applying the harmonic projection and leastsquares method, recovers most of the weaker Kelvin waves in the K a priori signal and reproduces their approximate timing and intensity (Fig. 5c).In the few locations where visible discrepancies between We now consider the error that is present in the reconstructed signal K reconst relative to the original signal K a priori , specifically = K reconst −K a priori .When the 100-year artificial timeseries is analyzed, it is found that the normalized root-mean-square (RMS) error 2 1/2 / K 2 reconst 1/2 is dependent on location along the waveguide as well as whether the fields are spatially-and temporally-filtered 255 (Figure 6a).(Here the angle brackets denote temporal averaging over the entire 100-year time span of the simulation, but no spatial averaging other than filters applied prior to the error calculation.)The error in recovering the original Kelvin wave signal is highest near the equatorial-coastal transition of the waves, and on the eastern end of the domain; elsewhere it is confined to a fairly narrow range.However, the error magnitude also depends on whether a spatial or temporal averaging filter is 260 applied prior to the error calculations.Except for the most error-prone regions, the error associated with unfiltered pointwise values of K is 50% to 60% of the total standard deviation of K.If the K has a spatial moving average (boxcar) filter applied, but temporal averaging is limited to 10 day ranges (the resolution of the original points), the normalized error decreases slightly in most locations and is smoother across the waveguide.The error associated with 30-day moving averages The probability and cumulative distribution functions associated with errors in K illustrate that errors of the same magnitude as the Kelvin waves themeselves are infrequent when a 30-day moving average filter is applied (Fig. 6b).Relative to the total standard deviation in filtered K, σ K , 270 the magnitude of the errors only exceed 0.5σ K about 10% of the time (either positive or negative), and only exceed 1σ K about 2.5% of the time.In this simulation, σ K ≈ 1.9 × 10 4 m 2 , so the error magnitude is less than 1×10 4 m 2 over 90% of the time.If the error is considered relative to the magnitude of the filtered reconstructed K reconst at each location and time, the error variance is somewhat larger.Even so, with the weakest Kelvin waves (|K reconst | < 0.3σ K ) excluded, the error will only 275 result in misdiagnosing the sign of most Kelvin waves (i.e., /K reconst > 1) approximately 8.5% of the time.Moreover, the 8.5% incidence decreases further if the threshold for excluding weak Kelvin waves is raised (approximately 5.5% for a 0.5σ K threshold and 2% for a 1σ K threshold); thus sign misdiagnosis using this method is rarely a problem for moderate and strong Kelvin waves.

Kelvin wave-related and unrelated SLA characteristics 280
For a more comprehensive view of the variability encapsulated by the Kelvin wave coefficient K described here, we also consider the reconstructed Kelvin wave-associated SLA field h K , in the context of the total SLA field h SLA observed by satellites over the 21-year period from 1992 to 2013.The SLA at each point in space and time can be considered as the sum of the reconstructed h K and a residual h res 285 with h res in theory encompassing contributions to the SLA field from processes unrelated to Kelvin waves, including Rossby and other planetary waves.Figures 7a,b illustrates the variances of h K and h res respectively, normalized by h SLA .The variance ratios of both h K (Fig. 7a) and h res (Fig. 7b) to h SLA exceed 1 along some parts of the waveguide, though this is more commonly the case with h res .
(NB: The variance ratios can exceed 1, since the h K and h res fields are generally not orthogonal.

290
Thus we do not describe the variance ratios as "explained variance" in the traditional sense, but rather compare the variances attributed to Kelvin waves vs. the residual to examine whether the residual signal is consistent with other phenomena such as Rossby waves.) Additionally, we compute the correlation between h K and the total h SLA and h res fields (Figure 8), to consider whether the sign of Kelvin wave activity covaries with that of the total SLA field and the waveguide (Fig. 8a) is strongly positive in the eastern part of the basin and along the coast, but insignificant or negative in the western part of the basin.The correlation of h K to h res is negative over nearly the entire domain, suggesting the tendency of h K and h res to be of opposite sign and explaining how h res in particular can have a much larger variance than the total SLA field.
The variance ratios (Fig. 7) and correlations (Fig. 8) suggest different contributions from h K and 305 h res variability in at least four distinct regions along the waveguide.In the western and central parts of the equatorial basin, even the maximum variances of h K near the equator are only slightly more than half the variance associated with h res .In this equatorial region, it is likely that most of h res can be attributed to Rossby waves; indeed a linear wind-forced model of equatorial waves (Nagura and McPhaden, 2010) has shown that in the western part of the basin, Rossby waves are associated with a 310 higher SLA standard deviation than Kelvin waves, even at the equator where Kelvin wave variability peaks.The correlation of h K and h res (Fig. 8) is also strongly negative here, consistent with the expected response of the ocean to a uniform zonal wind forcing, which would generate Kelvin and Rossby waves of opposite sign.In the eastern part of the basin and along the coast of Sumatra, the variances of h K and h res are more comparable, though this does not quite resemble the results from 315 the linear forced model (Nagura and McPhaden, 2010) which show a much larger component of SLA due to Kelvin waves than Rossby waves.Near the coast of Java, the variance of h K is much larger than that of h res , suggesting that most of the SLA variability in this area can be attributed to coastal Kelvin waves.Along Nusa Tenggara (NT) the variance of h res is once again comparable to or greater than h K ; this may be due in part to the complexity of the islands and bathymetry here.

320
It is also likely that the computation of K does not accurately resolve the splitting and diversion of Kelvin wave energy through Lombok Strait between Java and NT (e.g., Syamsudin et al., 2004;Drushka et al., 2010), since the least-squares fit exhibits a preference for slow tapering of K rather than the abrupt change in wave activity associated with the narrow strait.
Finally, the lack of a robust correlation between h K and h SLA along the equator in the western part 325 of the basin (Fig. 8) implies that using raw SLA to track Kelvin wave propagation may not accurately represent where waves originate.Namely, SLA crests and troughs that only become clearly apparent in the eastern part of the basin may actually have origins further west; some specific cases of this are evident when comparing the SLA and K values for 1997 (Fig. 3a, c).

Conclusions 330
The harmonic projection and least-squares method outlined here produces a measure of Kelvin wave activity that can be applied directly to satellite observations of SLA, not only along the equator, but also along low-latitude coastal waveguides.The method removes the westward-propagating signals tend to be of the opposite sign of the Kelvin waves on which they are superimposed.Therefore the use of this method helps to isolate the Kelvin wave-associated SLA signal; it also allows for some variation in the phase speed of the waves, so a comparison of K values with the results of linear wind-forced models of equatorial waves (e.g., Yu and McPhaden, 1999;Nagura and McPhaden, 2010) may be useful in studying some weakly nonlinear aspects of Kelvin waves.

345
It is notable that the SLA along the equator in the western Indian Ocean is not robustly correlated with Kelvin wave activity deduced from this method, a result that has important implications for the interpretation of SLA variability at the equator.The use of this method helps show that the SLA signal of numerous intraseasonal Kelvin waves can be traced to the western equatorial Indian Ocean and co-located with zonal wind stress forcing (e.g., Delman et al., submitted); therefore the lack of 350 readily identifiable eastward-propagating sea level anomalies at the equator at a given time does not necessarily imply that Kelvin waves are absent.Rossby waves may be obscuring the Kelvin waves' signal on the western side of the basin, and the computation of K may assist in tracking Kelvin waves from their true generation region.
We also observe two limitations of the harmonic projection/least-squares method in the form 355 presented here, and consider how these might be addressed.The first is the difficulty of resolving Kelvin wave activity to the east of Lombok Strait (∼115 • E), as evidenced by the abrupt increase in residual variance h 2 res (Fig. 7b) and decrease in h K -h SLA correlation (Fig. 8a) from the Java to the NT coastline.Prior analyses of altimetric SLA (Syamsudin et al., 2004;Drushka et al., 2010) indicate that about 30-50% of the Kelvin wave energy from the Java coastline continues eastward 360 along the NT coastline; the rest of the energy is presumed to transit north through Lombok Strait.In terms of our method, which does not track Kelvin wave energy through the strait, this would require an abrupt "dissipation" of Kelvin wave activity at Lombok Strait, which is likely not resolved by our tapered basis functions.Moreover, the overall skill of the method decreases approaching the eastern boundary of our domain and the complex topography of the Savu Sea region.One possible way to 365 resolve the abrupt change in Kelvin wave energy across the strait using our method is to adjust the weighting w of the misfit of the vector in (8).Namely, the misfit for the tapers that span the Lombok Strait region could be weighted more heavily, so that there is less tendency for the model to continue steady propagation of Kelvin waves past the strait.Additionally, a finer tapering scale (e.g., ∆x = 300 km instead of 600 km) could be adopted in this particular region, though errors due to altimetry Ocean Sci.Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -1, 2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.
The second issue is that this method was developed with the specific goal of isolating Kelvin waves, while the SLA signal due to Rossby waves is treated as a residual.By contrast, linear windforced models and prior studies that have used SLA projections have sought to quantify the activity of Kelvin and gravest-mode equatorial Rossby waves simultaneously.The variance ratios (Fig. 7) 375 and correlations (Fig. 8) suggest that the majority of the residual SLA field along the equator can be attributed to Rossby wave activity, and in principle there is no reason why our method could not be expanded to specifically target Rossby waves as well.One possible alteration to our method is to carry out the y-projections and the x,t-projections simultaneously; i.e., link each propagating basis function to the y-profile of a Kelvin or Rossby wave mode depending on its phase speed.Isolating 380 the SLA displacement associated with each mode would allow for a more complete picture of equatorial dynamics to be deduced from satellite altimetry, and this altered method will be explored in a future study.
145x s = x W the western boundary, while for s > 1, x s = x W + (s − 1)∆x.The forcing and dissipation of a wave within the domain can be expressed as the superposition of basis functions with varying s-values.OceanSci.Discuss., doi:10.5194/os-2016-1,2016   Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.Furthermore, to reduce the number of basis functions and make the subsequent least-squares problem less underdetermined, we limit the basis functions to certain phase speed ranges associated 150 with the waves we expect to observe using satellite altimetry.Therefore only basis functions A m,n,scorresponding to phase speeds c m,n = f n /k m typical of Kelvin waves (1.5 m s −1 ≤ c m,n ≤ 3.5 m s −1 ), Rossby waves (−1.2 m s −1 ≤ c m,n ≤ −0.4 m s −1 ), and stationary signals (k m = 0 or f n = 0) are included in A m,n,s , while the other basis functions are excluded.This phase-speed limitation reduces the number of basis functions to approximately twice the number of K y values in b.The 155 tapered basis functions in A m,n,s corresponding to the phase speed ranges are projected onto the vector b containing the K y values.First, a least-squares planar fit to b is removed, and then the basis functions are projected onto b m P = A T b (6) producing a vector m P of data projections that can be used to solve for the basis function coefficients m of identical size.For the projection values to be unbiased with respect to s, the projections are 160 normalized by the size of the nonzero domain, expressed as a normalizing vector n P , with n P m,n,s = (2/N )x max [x max − (x s − 0.5∆x)] −1 , and N = length(b).Accordingly, the normalized vector of basis function projections is m nP = n P m P = n P A T b (7) 2.4 Least-squares optimization and removal of westward-propagating signals After the x-t projections have been carried out and normalized, basis function coefficients are re-165 covered from the data projections in m nP using least-squares methods; the solution is optimized to prevent a poorly-scaled solution with the cancellation of large coefficients in m.In addition to constraining the size of the basis function coefficients m T m, we chose to minimize the misfit of the rate of change in data projection values along the waveguide, ∂m nP /∂s, in order to constrain high-wavenumber variability within the domain.We also minimize the misfit of the untapered data 170 projection values, m nP | s=1 , to the s = 1 basis function coefficients.Hence the vector that we minimize the misfit for is OceanSci.Discuss., doi:10.5194/os-2016-1,2016   Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.
e., sgn(f n ) = −sgn(k m ), are set to zero.The resulting vector m K is used to reconstruct the K y field with the westward-propagating signals removed b K = Am K (11) where vector b K consists of the Kelvin wave coefficients K as a function of x and t.These computations are carried out for overlapping two-year subsets of the data, which are then merged to create 185 a continuous field of K values for the period Sept. 1992-Dec.2013 covered by the alongtrack SLA dataset.3 Representations of Kelvin wave activity and error/variance estimates 3.1 Comparison of K values with raw SLA To demonstrate how well K represents Kelvin wave activity, we present a case study where we 190 compare the raw SLA along the IO equatorial-coastal waveguide during the year 1997 to the K y westward-propagating signals(e.g., Jan.-Feb., and Oct.-Dec.)unrelated to Kelvin waves, and most 200 likely represent Rossby waves flanking the equator.These westward-propagating signals are no Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.longer visible in the K values for 1997 (Fig. 3c), and the trajectories of alternating upwelling and downwelling Kelvin waves are much more readily apparent.The values of K can also be re-projected back into two spatial dimensions, to reconstruct the component of the SLA field that is associated with Kelvin wave activity.The reconstructed h K is 205 obtained by obtaining the wave amplitude h 0 associated with K 225ization (e.g.,Gentle, 1998) to construct m fields whose local covariance statistics in wavenumberfrequency (k-f ) space resemble values computed from the altimetry data, so that realistic Kelvin and Rossby wave signals could be generated.The m coefficients were adjusted so that their values are partially dependent on the local wave amplitude at the same wavenumber and frequencym| k,f,s = C s s−1 s =1 m| k,f,s + r(13) with C s the location-dependent adjustment parameter and r the Cholesky decomposition-based 230 random component.The variances of the basis functions were also adjusted so that the distributions Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License. of total Kelvin and Rossby wave variance along the waveguide are consistent with the variances obtained from satellite altimetry.Finally, after the artificially-generated eastward-and westwardpropagating signals were combined, a small amount of white noise was added to the K y fields; the variance of this noise is location-dependent and based on the variance in altimetry observations that 235 could not be explained by either Kelvin or Rossby wave signals.

K
a priori and K reconst are present (e.g., the intensities of the Kelvin waves in March-April, east of 90 • E), high amplitude westward-propagating signals and/or sharp noisy gradients are present in the 250 K y field.

295
the residual.The effective degrees of freedom N * for the correlations in Fig.8were computed from the decorrelation timescales of h K .With 21 years of data, values of N * range from approximately 50 to 500 over the spatial domain, with decorrelation timescales ranging from intraseasonal to semiannual.(For N * = 50, correlation coefficient magnitudes exceeding 0.23 exceed the 95% confidence threshold; for N * = 500 this threshold is approximately 0.08.)The correlation of h K to h SLA along 300 10 Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.
(i.e., Rossby wave signals) along such waveguides, and produces K coefficients that represent the time variability of Kelvin wave activity at each location along the waveguide.When filtered to re-335 11 Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.move sub-monthly temporal variability, the values of K have an RMS error of approximately 0.4 times the local standard deviation of K at most points along the waveguide; excluding the weakest waves, the method also diagnoses the sign of the Kelvin waves correctly over 90% of the time.A decomposition of the near-equatorial SLA into Kelvin wave-associated and residual components generates a residual field generally consistent with the activity of wind-forced Rossby waves, which 340 Acknowledgements.Andrew Delman (ASD) was supported by a NASA Earth and Space Science Fellowship, grant number NNX13AM93H.Janet Sprintall (JS) and Julie McClean (JLM) were supported by NASA 385 award number NNX13AO38G.Lynne Talley (LDT), JLM, and ASD were also supported by NSF grant OCE-0927650.The altimeter products were produced by Ssalto/Duacs and distributed by Aviso, with support from CNES (http://www.aviso.oceanobs.com/duacs/).Computations were carried out on the Geyser cluster within the Yellowstone computing environment hosted by the National Center for Atmospheric Research (NCAR).13 Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.Yuan, D., Rienecker, M. M., and Schopf, P. S.: Long-wave dynamics of the interannual variability in a numerical 430 hindcast of the equatorial Pacific Ocean circulation during the 1990s, J. Geophys.Res., 109, C05019, 2004.Ocean Sci.Discuss., doi:10.5194/os-2016-1,2016 Manuscript under review for journal Ocean Sci.Published: 12 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 1 .
Figure 1.AVISO gridded sea level anomaly (SLA) in the Indian Ocean basin, on (a) June 25, 1997 and (b) July 16, 1997, with upwelling Kelvin waves (depressed SLA) peaking in the central Indian Ocean and along the Sumatra/Java coasts respectively (green dashed ellipses).The brown dashed lines indicate the radii of deformation for 1st baroclinic mode Kelvin waves, with the radius extended along the Indonesian coastal waveguide south of the equator.The locations of Sumatra, Java, and Nusa Tenggara (NT) are also indicated.

Figure 2 .
Figure 2. Schematic illustrating the use of tapered basis functions.(a) Profile of a non-tapered two-dimensional basis function Am,n in x and t.(b) Profile of a tapered basis function Am,n,s, with tapering location x = xs and a tapering window of ∆x.(c) Profile of two superimposed basis functions Am,n,s − A m,n,s , with tapering locations of x = xs and x = x s respectively; the superposition of two or more tapered basis functions allows for the forcing and dissipation of waves.

Figure 3 .
Figure 3. (a) Sea level anomaly (SLA) along the waveguide that crosses the equatorial Indian Ocean and follows the coasts of Sumatra, Java, and Nusa Tenggara, during 1997.The data plotted are from the 1/3 • gridded AVISO product, in units of cm.(b) Kelvin wave y-projections Ky along the waveguide during 1997, derived from the alongtrack SLA data collected by the TOPEX/Poseidon satellite during 1997, in units of 10 4 m 2 .(c) The Kelvin wave coefficient K values during 1997, in units of 10 4 m 2 .The vertical dashed lines in each plot indicate (from left to right) the locations of the equatorial-coastal transition (98 • E), the Sumatra-Java transition at Sunda Strait (105 • E), and the Java-Nusa Tenggara transition at Lombok Strait (115.65 • E).18

Figure 4 .
Figure 4. Maps of (a)-(e) SLA and (f)-(j) reconstructed hK , the Kelvin-wave associated SLA, for snapshots (dates in the top-right corner of each panel) over a 2-month period in 1997.As in Fig. 1, the brown dashed lines indicate the radii of deformation for 1st baroclinic mode Kelvin waves.

Figure 5 .
Figure 5. (a) Kelvin wave y-projection values Ky from year 5 of the Monte Carlo simulation, along the Indian Ocean equatorial-coastal waveguide.(b) Kelvin wave coefficient values Ka priori generated in year 5 of the Monte Carlo simulation.(c) Reconstructed Kelvin wave coefficient values Kreconst for year 5 of the Monte Carlo simulation, obtained using the harmonic projection and least-squares method.All quantities are given in units of 10 4 m 2 .

Figure 6 .
Figure 6.(a) Normalized root-mean-square (RMS) error estimates for K as a function of longitude along the waveguide, based on the 100-year Monte Carlo simulation and computed as 2 1/2 / K 2 1/2 .The different curves show the effect on RMS error of applying spatial and temporal moving average filters to the a priori and reconstructed K values.(b) The probability distribution function (solid lines) and cumulative distribution function (dashed-dotted lines) of normalized error, with a 20 • longitude and 30 day moving average filter applied prior to the error computation.The curves are shown for two different normalizations: (blue) normalized by the standard deviation of K over all longitudes and times /σK , and (green) normalized by the filtered reconstructed value Kreconst for each longitude and point in time /Kreconst; in the latter case errors associated withe values of |Kreconst| < 0.3σK have been excluded.

Figure 7 .
Figure 7. (a) The variance ratio of Kelvin wave-associated SLA to total SLA, h 2 K / h 2 SLA .(b) The variance ratio of the residual to total SLA, h 2 res / h 2 SLA .The annual and semiannual harmonics have been removed by linear regression.

Figure 8 .
Figure 8.(a) The correlation coefficient of Kelvin wave-associated SLA to total SLA, hK hSLA / h 2 K h 2 SLA 1/2 .(b) The correlation coefficient of Kelvin wave-associated SLA to the residual, hK hres / h 2 K h 2 res 1/2 .The annual and semiannual harmonics have been removed by linear regression.Only correlation coefficients exceeding the 95% confidence threshold for significance (based on a Student's t-distribution) are shaded.