Harmonic analysis of tide gauge data in the north-eastern Bay of Bengal: Amplitudes and residuals for 366 day intervals

Figure 1 shows a map of the head bay of the Bay of Bengal and the location of 8 tide gauges. Here we refer to this region simply as north-eastern Bay of Bengal. UHSLC maintains an archive of research-quality water level data with hourly and daily resolution, containing (among many others) the depicted gauges. Note that the locations of stations "Khal #10" and Chattogram (formerly known as Chittagong) are barely distinguishable on the map. Khal #10 seems to be located slightly upriver of the Chittagong station (detail not shown). Note also that the location coordinate metadata for the tide gauge in Cox's Bazaar may be incorrect. A detailed description of the bathymetry around individual stations will be presented in a future study along with results from numerical modelling.

The tide gauge data is preprocessed by UHSLC (Caldwell et al., 2015) and/or possibly other third parties. IOC (1992) describe the standardized data processing and quality assurance methodology followed to produce the Research Quality dataset archived by the NOAA National Centers for Environmental Information as the Joint Archive for Sea Level (JASL, see Caldwell et al. 2015). The quality assessment policy (FTP link) states that assessment is mostly based on residuals (observed data minus predicted tides) of hourly data. IOC (1992) state that the hourly data is computed from the high-frequency data delivered by the instrument via a three-point Hanning filter. In at least some of the station metadata described below, it is explicitly mentioned that the filter is centered on the hour. IOC (1992) further state that the filter was chosen to minimize aliasing from subsampling and to minimize the attenuation of the tides. The archive also contains daily time series. Caldwell (2015) describes the two-step tidal filtering procedure for computation of the daily values. In a first step, dominant diurnal and semi-diurnal tidal components are removed from the quality controlled hourly values. Secondly, a 119-point convolution filter is applied to remove the remaining high frequency energy and to prevent aliasing. In this post we mostly use hourly data, but also show filtered ("de-tided") daily data in the plots of the monthly mean sea levels and the residual time series.

In this blog post we present results from harmonic analysis and discrete Fourier transforms, and point out that we do not assess the sensitivity of our results to the choice of any low-pass filters used to produce the hourly data. In particular, any energy contained in high frequencies of the physical (continuous) quantities, which may not have been filtered, will project onto discretized lower frequencies via aliasing.

Each station (i.e. time series) is accompanied by an individual quality assessment report issued by UHSLC, containing information such as missing data, notes on the reference level, spurious fluctuations, etc. (see e.g. here for examples for the Indian Ocean, but note that this is an FTP link which your web browser may refuse to follow ). Each report includes

• a Completeness Index ("CI") indicating the percentage of hours with available data, which may include replacements of gaps of less than 25 hours by prediction data
• a list of gaps longer than a day. These are not patched with prediction data, and instead are filled with a missing data flag. The format of the list is e.g. "44-45,325-330,352-353" referring to 3 gaps, the first of which starts at some (unspecified) hour on year-day number 44 and ends at some hour on year-day number 45. We found that the exact hours can most easily be retrieved by searching for the missing data flag in a human-readable format of the data files (e.g. the comma-separated-value format distributed by UHSLC).
• a summary of days with consecutive missing or obviously wrong data in a time span longer than 6, but less than 25 hours that have been replaced with prediction data, with hours in parentheses.
• a summary of days with questionable fluctuations in the residual series. Fluctuations within the residuals are considered significant and are noted if they are greater than 25 cm. This policy has changed slightly as of October 1992: only suspicious fluctuations due to obvious instrumentation problems are noted. Minor timing drifts are no longer included in the list of questionable days.

For our purposes we pick out years with a CI of 100% and with as few questionable fluctuations as possible. We do not consider the stations Khal #10 and Teknaf. Some time series contain replaced gaps. This implies that our harmonic analysis is based on a previous harmonic analysis performed by UHSLC, and hence could be described as somewhat circular or iterative. One motivation to repeat the analysis is simply that UHSLC do not seem to publish their computed harmonic constants, perhaps because they are merely a by-product of their quality control procedure. Also note that one particular advantage of harmonic analysis based on least squares over discrete Fourier analysis is the ability to handle data gaps and irregularly spaced data. Since there are no data gaps in our (already patched) time series, and the interval is constant, the question arises what further advantages harmonic analysis has. This is issue not further discussed in this blog post and the reader is referred to the literature (e.g. Parker, 2007, their section 3.4.5).

The general method of harmonic analysis based on the least squares solution is described e.g. in Foreman and Henry (1989).

Given an observed time series \(\zeta(t)\) of water level above a reference datum, one may define a predicted water level as

where \(n\) is the number of tidal constituents and \(\omega_i\) the constituent's angular velocity. The objective is to find harmonic constants, i.e. amplitude \(a_i\) and phase \(\phi_i\) for each constituent (and the constant \(Z_0\)), such that the residual time series

is minimal in a least-square sense. Amplitudes and phases may be modified by nodal corrections.

We use the SLP64 software package distributed by the University of Hawaii Sea Level Center, which includes the IOS Tidal Package written by Foreman (1977). The IOS Tidal Package in turn is based on the work by Godin (1972) and Godin and Taylor (1973). The related T_Tide MATLAB package of Pawlowicz et al. (2002) is based on the IOS Tidal Package but includes methods to compute confidence intervals and signal-to-noise ratios. Originally, a secondary motivation of this blog post was to collect useful notes on the use of the IOS Tidal Package and T_Tide and conflate them with a practical application to data from the Bay of Bengal. However, the blog post became too long and we moved the software notes into a separate post.

The notes on IOS Tidal Package and T_Tide describe some differences between the software packages in more detail. Here we only point out that there are slight differences in the predicted water levels, and hence the residual time-series. First, the packages apply nodal modulation corrections differently in their tidal predictions. Secondly, T_Tide is able to compute confidence intervals or signal-to-noise ratios for each harmonic constant (Pawlowicz et al., 2002). In its default configuration, T_Tide version v1.3_beta will include only significant constituents in the prediction of the water level, and hence also in the computation of the residual water level. More details are given in the companion blog post.

For this blog post we simply adopt the default parameters of the t_tide.m routine in T_Tide version v1.3_beta. At this stage we do not attempt to justify this choice in the context of this study. Instead, the results presented here serve to gain intuition about the meaning and usefulness of these confidence intervals. Hence the figures of residual time series below will show two residuals throughout this blog post, where one is computed with IOS Tidal Package using all constituents, and the other is computed with T_Tide using only significant constituents. We interpret a constituent as statistically significant if its harmonic constants have an SNR>2, following the default configuration of T_Tide.

In applying the harmonic analysis method to the UHSLC data we loosely follow the steps outlined in chapter 4.5 of Parker (2007) and the manual of SLP64. We analyze time series of 366 consecutive days, as outlined in the manual of SLP64. For non-leap years, we choose the 365 days of the calendar year, and add the first day of the subsequent year. According to the constituent sorting strategy of Foreman (1977), the number of resolvable constituents with this length is 67, i.e. 44 main constituents (excluding the constant constituent Z0) and 24 shallow-water constituents. Note that the astronomical constituent Γ2 listed in Table 1 of the companion blog post is not resolvable, since a time series of about 472 days would be necessary to separate it from H1 (α2) (Foreman, 1977, their Table 3)..

As an example for the data processing, consider the year 2008. This is a leap year and hourly observations from 2008/1/1 00:00 UTC to 2008/31/12 23:00 UTC yield 366*24=8784 observations. Both IOS Tidal Package (provided with SLP64) and T_Tide will clip the last observation to obtain an odd number of observations (Foreman, 1977, their section 2.2.1).

Amplitudes of tidal constituents:

• Chittagong 124a
• Hiron Point 134a
• Cox's Bazaar 136a
• Charchanga 138a
• Khepupara 139a
• Akyab (Sittwe) 907a

Monthly mean sea level and predictions by the solar annual constituent:

• Chittagong 124a
• Hiron Point 134a
• Cox's Bazaar 136a
• Charchanga 138a
• Khepupara 139a
• Akyab (Sittwe) 907a

Residuals of harmonic analysis:

Figure 2 shows amplitudes of tidal constituents for multiple years at Chittagong (124a). The (K1+O1)/M2 mean amplitude ratio is 0.16 and the tide is generally of semidiurnal type (c.f. Parker, 2007, their page 43). At all other stations this quantity does not exceed 0.25. M2 at Chittagong (1.73m) is much larger than at the other stations. M2 ranges from 0.76m (Akyab) to 0.976 (Charchanga) at all other stations. Regarding the seasonal cycle of sea level in the region, note the relatively large amplitude of SA, whose multi-year average is more than 0.3 m. At the other stations the multi-year average of SA is 0.345 m (Hiron Point), 0.398 m (Cox's Bazaar), 0.485 m (Charchanga), 0.398 m (Khepupara) and 0.268 m (Akyab). In 4 out of 6 stations (i.e. Hiron Point, Cox's Bazaar, Charchanga and Khepupara), the mean SA exceeds the mean S2, and thus SA has the second largest mean amplitude of all constituents. Contrary to M2 and other constituents, SA is not much larger at Chittagong. This may be expected for any long-period constituent, which vary in general on larger spatial scales. The constituents H1 and H2 are 0.150 and 0.099, respectively. In Section 5 we discuss the relation of H1 and H2 to the constituents \(MA_2\) and \(Ma_2\) introduced by Corkan (1934) to represent annual modulations of M2 due to meteorological phenomena (see also Müller et al., 2014).

Figure 3 shows monthly mean sea level anomalies for various years at the Chittagong station. Plots for the remaining stations can be accessed from the list. The multi-year average has a range of about 0.7 m, and is reasonably well approximated by the sinusoidal contribution of the solar annual constituent SA, with an amplitude of 0.339 m (c.f. Figure 2) and an argument computed from the multi-year average Greenwich epoch and an arbitrary equilibrium argument (here for January 1, 2008). Figure 18 shows the monthly means for the station at Charchanga, where the average seasonal range is more than 0.8 m (c.f. Fig. 2 of Tazkia et al., 2017).

Figure 4 shows \(\zeta_{res}\) at the Chittagong (124a) station for the year 2008. All other stations and years can be accessed via the list above. As noted in Section 3 , we choose to show two different residuals throughout this blog post, one calculated with IOS Tidal Package and T_Tide. This is to highlight any possible differences arising from (1) the exclusion of non-significant constituents (as done by T_Tide in its default configuration), and (2) the different application of nodal correction numbers for predictions, as described in the companion blog post. The black line in Figure 4 shows residual heights at Chittagong for 2008, calculated with the IOS Tidal Package. The dashed red line shows T_Tide predictions omitting non-significant constituents. The dashed green line shows filtered daily values obtained directly from UHSLC.

Apart from transient fluctuations on various temporal scales, the residual obviously oscillates with frequencies in the semidiurnal band. The variance of the residual is 0.023 \(m^2\). Table 12 lists the multi-year mean of the residual variance for all stations.

The variances of the observed, predicted and residual time series for the year 2008 at Chittagong are

• observed: 1.840 \(m^2\)
• predicted:
  • IOS Tidal Package: 1.819 \(m^2\)
  • T_Tide: 1.813 \(m^2\)
• residual:
  • IOS Tidal Package: 0.023 \(m^2\)
  • T_Tide: 0.027 \(m^2\)

The variances of the observed and predicted time series are shown as function of the tidal species in Figure 5. The line for \(\zeta_{pred}\) is computed with IOS Tidal Package (all constituents). The units of the horizontal axis are frequencies corresponding to the period of a mean lunar day, i.e. about 24.8412 hours (c.f. Parker, 2007, their table 2.1. on p.29). The variance shown is accumulated within each band. The observed variance has a maximum of 1.742 \(m^2\) at the semidiurnal band, and is about 0.03 \(m^2\) in the diurnal band. It is noteworthy that the variance is generally higher at even multiples of \(f_L=\omega_L/(2\pi)\). We speculate that this is because there is far more energy in the semidiurnal band than in the diurnal band (both in the water level and the residual, as shown below), and that the energy at higher frequencies is dominated by overtides of the semidiurnal species.

Figure 6 shows how the variance of the residual is distributed over the tidal species. The accumulated variance in the low-frequencies (up to 0.5*\(\omega_L\)) is about 0.012 \(m^2\), and about 0.008 \(m^2\) in the semidiurnal band. Like for the observed and predicted time series, the variance of the residual is generally higher for even multiples of the lunar frequency. However, the energy in the residual is highest in the low-frequency band.

Figure 7 shows a discrete Fourier transform of the residual time-series in the low-frequency band. There is a lot of energy between the semi-annual and monthly frequencies, presumably due to meteorological phenomena. Note that the resolution in frequency space is about 1/y. Higher resolution of the discretization would require time series of larger intervals than one year.

Figure 8 shows a discrete Fourier transform of the residual in the semi-diurnal band. There is some energy in frequencies which could be resolved with time-series of 366 days. For example, the shallow-water constituent ST2 is in between the resolved constituents \(\epsilon_2\) (with shallow-water equivalent \(MNS_2\), see table A.1. of Parker 2007) and the lunar constituent \(2N_2\) (shallow-water equivalent \(2NM_2\)). The frequency interval \(ST2-\epsilon_2\) and (\(2N_2-ST2\)) is larger than required by the Rayleigh criterion, as can be seen from Figure 8 by comparison to the interval between any two narrowly spaced black vertical lines indicating resolved frequency pairs. Configuring T_Tide to include 2NS2, ST2 and 2SM2 decreases the total residual variance from 0.023 \(m^2\) to 0.021 \(m^2\), and has no clear visible effect on the residual time series (not shown).

In the context of the second method used by Müller et al. (2014), i.e. the method involving H1 and H2, two questions arose during the writing of this blog post. The first relates to the exact definition of H1 and H2, the second arises from what appears to be a mistake in Eq. (22) of Müller et al. (2014). The first issue is discussed in the remainder of this section. The second issue is discussed in Appendix B.

We noted an apparent inconsistency in the naming conventions of tidal constituents. The solar annual constituent named SA seems to be defined alternately as the inverse of a tropical year and an anomalistic year. For example, in table A.1 of Parker (2007), SA has Cartwright numbers (0,0,1,0,0,0) and hence the period of a tropical year (see also our table 23 ). In general, governmental agencies in the U.S. seem to use the tropical year for SA. As an example, see the listing of tidal constituents on the web page of the Center for Operational Oceanographic Products and Services (CO-OPS). On the other hand, in the IOS Tidal Package and T_Tide, SA has Cartwright numbers (0,0,1,0,0,-1) and a frequency of 0.0001140741 cycles per hour, which represents an anomalistic year (to machine precision). Furthermore, in table 1(a) of Cartwright and Edden (1973), the Cartwright number (0,0,1,0,0,0) does not appear, the closest frequency vectors being (0,0,1,0,0,-1) and (0,0,1,0,0,1). This inconsistency seems to carry over to the "meteorological constituents" proposed by Corkan (1934), who stated that the "arguments" of M2 on the one hand, and his newly introduced \(MA_2\) and \(Ma_2\) on the other hand, differ by a magnitude of \(h\), which denotes the mean longitude of the sun. E.g. Schureman (1958) defines the term tropical year in their table 1 (on their p. 163) as

Their table 1 contains a note stating that (quote:) "symbols refer to the rate of change in mean longitude", and further above in their table 1, the symbol \(h\) is associated with the sun. It seems that the usage of \(h\) in both Schureman (1958) and Corkan (1934) is consistent, and hence that Corkan (1934) did not suggest the involvement of the solar perigee. This is inconsistent with the naming convention of Müller et al. (2014), who state that the annual constituents proposed by Corkan (1934) are sometimes referred to as \(\alpha_2\) and \(\beta_2\) in the literature. We were unable to obtain their reference ( Hartmann and Wenzel 1995) free of charge, but e.g. Parker (2007) includes at least \(\alpha_2\) in their tables A.1 and A.2, with Cartwright numbers (2,0,-1,0,0,1), i.e. in terms of frequency

where (from Parker, 2007, their table 2.1)

The inclusion of \(f_5\) seems to deviate from what Corkan (1934) proposed. Searching for further clarification, we found a tidal constituent list by the IHO distinguishing MA2 from alpha2 in the expected way, i.e. MA2 has Cartwright numbers (2,0,-1,0,0,0) and MB2/MA2* have Cartwright numbers (2,0,1,0,0,0), even though the latter are listed with then name M(KS)2, which is apparently a compound constituent. The document contains a note stating that

Despite their appearance neither [MA2 nor MB2], nor constituents of other species which include A and B, are compound constituents – there are no constituents A or B to form a compound. They are constituents whose speeds differ by one cycle a year from that of M2. The A in MA2 was intended to signify the Annual differences. MB2 was originally called Ma2 but this became ambiguous when spoken, or typed on computers without lower case, and so it was initially changed to MA2*. However, this in turn was thought to be clumsy and hence MB2 was finally adopted. Although theoretically they should have the same values of u and f as M2, they are so small that they are commonly treated as having values of u = 0 and f = 1.

Only after reading this paragraph, the author noted that the constituent Sa in the IHO list also has two definitions, i.e. (0,0,1,0,0,-1) and (0,0,1,0,0, 0). In other words, the apparent inconsistency in the constituent naming of Müller et al. (2014) may be a logical consequence of the inconsistency in the naming of SA.

Do the differences between MA2 and MB2 on the one hand, and \(\alpha_2\) and \(\beta_2\) (or its shallow water equivalent) on the other hand, matter in practical applications with time series of 366 days? Given the interval which would be necessary to separate the pairs, i.e. the period of the rotation of the perihelion around the sun (~20940 years), this may not be the case.