See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/249605162

Versatile Harmonic Tidal Analysis: Improvements and Applications

Article  in  Journal of Atmospheric and Oceanic Technology · April 2009

DOI: 10.1175/2008JTECHO615.1

CITATIONS READS

67 458

3 authors, including:

M. G. G. Foreman Josef Y. Cherniawsky

Fisheries and Oceans Canada Fisheries and Oceans Canada
134 PUBLICATIONS   6,459 CITATIONS    76 PUBLICATIONS   1,216 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

coastal aquaculture modelling View project

climate change View project

All content following this page was uploaded by M. G. G. Foreman on 25 January 2014.

The user has requested enhancement of the downloaded file.

806 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

Versatile Harmonic Tidal Analysis: Improvements and Applications

M. G. G. FOREMAN AND J. Y. CHERNIAWSKY

Institute of Ocean Sciences, Fisheries and Oceans Canada, Sidney, British Columbia, Canada

V. A. BALLANTYNE
Canadian Hydrographic Service, Fisheries and Oceans Canada, Sidney, British Columbia, Canada

(Manuscript received 13 December 2007, in final form 21 August 2008)

ABSTRACT

New computer software that permits more versatility in the harmonic analysis of tidal time series is
described and tested. Specific improvements to traditional methods include the analysis of randomly sampled
and/or multiyear data; more accurate nodal correction, inference, and astronomical argument adjustments
through direct incorporation in the least squares matrix; multiconstituent inferences from a single reference
constituent; correlation matrices and error estimates that facilitate decisions on the selection of constituents
for the analysis; and a single program that analyzes one- or two-dimensional time series. This new metho-
dology is evaluated through comparisons with results from old techniques and then applied to two problems
that could not have been accurately solved with older software. They are (i) the analysis of ocean station
temperature time series spanning 25 yr, and (ii) the analysis of satellite altimetry from a ground track whose
proximity to land has led to significant data dropout. This new software is free as part of the Institute of
Ocean Sciences (IOS) Tidal Package and can be downloaded, along with sample input data and an ex-
planatory readme file.

1. Introduction many of these codes are restrictive in both the form of

the input time series (e.g., regularly sampled, albeit with
There have been many advances in tidal analysis and
gaps) and the manner in which nodal correction, as-
prediction since the earliest documented predictions for
tronomical argument, and inference calculations are
the bore on the Chhien-Thang River in China and flood
made (e.g., as adjustments to results from a least squares
tide at London Bridge in the eleventh and thirteenth
fit). In the early days of harmonic analysis, these re-
centuries, respectively (Cartwright 1999). Parker (2007)
strictions arose from computer limitations that neces-
has recently published a guide on the various consid-
sitated efficiency more than accuracy in the algorithms
erations that are needed, and contemporary approaches
(Godin 1972; Foreman 1977, henceforth F77). However,
that can be used, to carry out accurate tidal analyses and
present computer capacities mean that these restrictions
predictions. One of the most successful and widely used
need no longer apply.
approaches has been, and continues to be, harmonic
In this study, we develop and test a more versatile
analysis wherein the energy at specific tidal frequencies
harmonic analysis technique that can accept randomly
is determined by a mathematical fitting procedure,
sampled data and embed the nodal and astronomical
usually least squares. Though computer software that
argument corrections and multiple inference calcula-
performs harmonic tidal analysis of one- and two-
tions into an overdetermined matrix that is solved using
dimensional time series has been are available for more
singular value decomposition (SVD) techniques (Golub
than 40 yr (links to software packages are available online
and Van Loan 1983; Press et al. 1992). The input time
at http://www.pol.ac.uk/psmsl/training/analysis.html),
series is also allowed to be one or two dimensions,
thereby eliminating the need for separate programs to
Corresponding author address: M. G. G. Foreman, Institute of
Ocean Sciences, Fisheries and Oceans Canada, P.O. Box 6000,
analyze tidal heights and currents. In the latter case, the
Sidney, BC V8L 4B2, Canada. final harmonics are also expressed in terms of current
E-mail: mike.foreman@dfo-mpo.gc.ca ellipse parameters. Though the use of randomly sampled

DOI: 10.1175/2008JTECHO615.1

Ó 2009 American Meteorological Society

APRIL 2009 FOREMAN ET AL. 807

data raises issues of constituent selection and indepen- solved this dilemma by defining constituent ‘‘clusters’’
dence, the SVD approach allows for the calculation of that have the same first three Doodson numbers. Each
covariance matrices and correlation coefficients that cluster was assigned the name of its major constituent
permit an assessment of these dependencies (Cher- (in terms of tidal potential amplitude), while the lesser
niawsky et al. 2001). Thus, an iterative approach can be constituents were termed ‘‘satellites.’’ Harmonic tidal
used to determine which constituents should be sought analysis then followed two steps: (i) all satellites were
directly and which should be inferred. Embedding the ignored and amplitude and phases were determined for
inference calculations into the overdetermined matrix all major constituents that could be resolved given the
means that the inferred constituents will affect all other time series length; (ii) a so-called nodal correction was
constituents included in the analysis, not only the ref- performed to account for the presence of the satellites
erence constituent that is the basis of the inference. In and—if necessary—an inference was carried out to
addition, embedding the nodal corrections in the matrix correct for important missing major constituents. More
has not only a similar effect for the satellite and major details on both these steps will be given later.
constituent, but it also removes the need for assump- The five astronomical variables referred to previously
tions that underlie usual postfit corrections and that may are associated with the 27.32 day, 365.24 day, 8.85 yr,
restrict the length of the analysis period (F77). 18.6 yr, and 21 000 yr cycles arising from variations in
The analysis is first tested by comparing its results the mean longitude of the moon, the mean longitude of
against those arising from the F77 conventional approach the sun, the longitude of the lunar perigee, the longitude
for a pair of synthetic time series with and without back- of moon’s nodal progression (inclination of its orbit to
ground noise. A direct assessment of accuracy is possi- the equator), and the longitude of the solar perigee,
ble, because the amplitudes and phases of the constit- respectively. Though the term ‘‘nodal correction’’ was
uents are known. A second test is then performed with a originally coined before the advent of modern com-
12-yr time series to demonstrate the effects of embed- puters to designate corrections for only the moon’s
ding the astronomical argument and nodal correction nodal progression that were not incorporated into the
calculations directly into the least squares fit, rather astronomical argument calculation for the main con-
than as postfit adjustments. To demonstrate the use of stituent, the term ‘‘satellite modulation’’ is more ap-
correlation coefficients in constituent selection, further propriate now, because the correction has been ex-
tests are performed with a synthetic time series of ran- tended to include the effects of variations in lunar and
domly sampled data. Finally, two examples are given to solar perigees. In mathematical terms, the harmonic
illustrate the versatility of the new technique. The first is analysis approach originally proposed by Godin (1972)
the analyses of salinities from conductivity–temperature– and employed by F77 and Pawlowicz et al. (2002,
depth (CTD) stations that span 20 yr, and the second is henceforth PBL02) assumed that a one-dimensional
the analysis of Ocean Topography Experiment (TOPEX)/ time series with tidal and nontidal energies can be ex-
Poseidon (T/P) satellite altimetry along a short track pressed as
crossing the Strait of Georgia, whose proximity to land
n
has led to significant data dropout.
h(tj ) 5 Z0 1 å f k (t0 )Ak cos [vk(tj  t0 )
k51
2. Traditional and versatile harmonic analyses 1 V k (t0 ) 1 uk (t0 )  gk ] 1 R(tj ), (1)
Tidal potential theory (Doodson 1921) predicts the
existence of hundreds of tidal frequencies, each of where h(tj) is the measurement at time tj; Z0 is a con-
which can be expressed as a linear combination of the stant background value; fk(t0) and uk(t0) are the nodal
rates of change of mean lunar time and five astronom- corrections to amplitude and phase, respectively, at
ical variables that uniquely specify the position of the some reference time t0 for major constituent k with
sun and moon. For each frequency, the six integer co- frequency vk; Ak and gk (k 5 1,n), are the amplitude
efficients associated with this linear combination are and phase lag of constituent k, respectively; Vk(t0) is the
referred to as its Doodson numbers. However, it is astronomical argument for constituent k at time t0; R(tj)
neither practical nor mathematically feasible to include is the nontidal residual; and n is the number of tidal
all constituents in every analysis, because many fre- constituents. A least squares approach is usually em-
quencies are so close that a time series of several years’ ployed to solve for Z0, Ak, and gk; the observation times
duration is required to separate some neighbors by one are often assumed to arise from regular sampling (e.g.,
cycle; while, according to potential theory, others hourly, though gaps are permitted); and the number n
should have very small amplitudes. Godin (1972) re- and specific constituents k selected for the analysis are
808 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

usually determined in accordance with the time series sumption is often invalid as a result of nonlinear in-
length and the estimated background noise level. teractions between the tide and storm surges or vari-
Deciding which major constituents should be in- able river discharge that literally changes tidal ampli-
cluded in the first stage of a harmonic analysis is not tudes and phases for the periods during which these
easy. Whereas the so-called Rayleigh criterion (Godin phenomena occur. However, this is generally not a se-
1972) argues that a time series of length T is required to rious problem, as the effects are often small or of short
distinguish between constituents with a frequency sep- duration. In cases where the effects are more substan-
aration of T 21, linear algebra suggests that m inde- tial, such as for internal tidal currents that change with
pendent observations h(tj) should be sufficient to solve a the stratification or seasonally varying ice cover that
matrix equation for m/2 amplitudes and phases, as de- can modify both tidal elevation and current harmonics,
scribed in Eq. (1). In actuality, the decision also needs to wavelet analysis (Jay and Flinchem 1997, 1999) has been
consider the nontidal signal [R(tj) in Eq. (1)] and the used successfully, though it does not actually produce
condition number (Ortega 1972) of the least squares improved amplitudes and phases.
matrix. Munk and Hasselmann (1964) refined the As described earlier, the second limitation is that the
Rayleigh criterion by showing that ‘‘meaningful state- time series may not be sufficiently long to adequately
ments’’ can be made about the tidal energies associated separate the energy from constituents that are close in
with frequencies v1 and v2 provided frequency, or are aliased as a result of infrequent sam-
pling. Though this can be overcome through the use of
jv1  v2 j . T 1 (signal/noise)1/2 , (2) inference—wherein the amplitude and phase of lesser
constituents are assumed to have specific relationships
where signal/noise is the ratio of the tidal variance to the to a larger reference constituent—the choice of infer-
nontidal variance. Foreman and Henry (1989) extended ence parameters and the manner in which the calcula-
the selection analysis further by using standard matrix tion is carried out can limit the accuracy. This will be
theory. Assuming that Ax 5 b and Ax9 5 b9 are the discussed and illustrated later.
matrix equations associated with (1) when the back- The third limitation arises from the implementation
ground noise R is zero and nonzero, respectively; K(A) of the nodal corrections and, to a much lesser extent, the
is the condition number for A; and k  k is a measure- astronomical argument. The use of vk(tj 2 t0) 1 Vk(t0)
ment norm, matrix theory (Ortega 1972) states that in Eq. (1) assumes that the astronomical argument for
constituent k, Vk(tj), varies linearly about some refer-
k x  x9 k k b  b9 k ence time t0. Though this is generally a very good as-
# K(A) . (3)
kxk kbk sumption, it does tend to break down as the time series
extends over several years. On the other hand, assuming
In the context of tidal analysis, this inequality has the that fk and uk remain constant with their values for t0 is
following interpretation: the term kb - b9k/k b k is the more questionable. For many constituents, not only do
inverse of Munk and Hasselmann’s (1964) signal-to- these nodal variations change significantly over 18.6 and
noise ratio, while kx - x9k/kxk is the fractional error in 8.85 yr but also the manner in which the corrections are
the fitted amplitudes and phases. The effect of including implemented (Godin 1972; F77; PBL02) can cause the
relatively close frequencies (in the Rayleigh criterion accuracy of analysis results to deteriorate as the time
sense) in the harmonic analysis is to make the rows of A series extend beyond one year.
more linearly dependent and increase K(A). So, the In an effort to remove some of these deficiencies and
combination of relatively close frequencies with sub- limitations, a new harmonic analysis program has been
stantial background noise relative to the signal should developed that starts by replacing Eq. (1) with
cause relatively large differences between the calculated n
set of parameters x9 and their true values x. However, if
the frequencies are not close and/or the background noise
h(tj ) 5 Z0 1 atj 1 å f k (tj )Ak cos [V k (tj ) 1 uk(tj)  gk]
k51
is small, the fitted solution should be more accurate. 1 R(tj ). (4)
Though Godin’s (1972) harmonic tidal analysis ap-
proach has been used successfully for many years, it In this case, V, u, and f are evaluated at the precise times
does have limitations. The first and perhaps foremost of each measurement, thus eliminating inaccuracies that
among these is the underlying assumption of statio- arise from assuming a linear variation in the astronomical
narity, wherein the tidal amplitude and phase for each argument [i.e., Vk(tj) 5 vk(tj 2 t0) 1 Vk(t0)] and tempo-
constituent are assumed to remain constant over the rally constant values for the nodal corrections. In addi-
period of the time series. In shallow water, this as- tion, this program includes a linear trend (coefficient a),
APRIL 2009 FOREMAN ET AL. 809

being accurate and efficient, SVD has other advantages

that will be explained shortly.
Though it might be viewed as a disadvantage rather
than an improvement, another change with this new
program is that the constituents to be used in the anal-
ysis are not selected automatically. They must be
specified by the user. This was deemed necessary, as the
provision of arbitrary sampling generally means that
variations of the Rayleigh criterion are no longer valid
for determining constituent selection. F77 and PBL02,
FIG. 1. Schematic showing the old and new harmonic tidal analysis for example, employ a Rayleigh criterion decision tree
approaches. (see Tables 1–4 in F77) that provides a hierarchy of
constituent selection based on frequency separation and
allows for the measurements tj to arise from arbitrary tidal potential amplitudes such that when a time series is
sampling, and permits multiconstituent inferences (i.e., not sufficiently long to separate two neighboring con-
more than one constituent can be inferred from a single stituents, only the one with the larger expected ampli-
reference constituent) that are computed directly within tude is included directly in the analysis. On the other
the least squares fit, rather than as a correction to postfit hand, the SVD approach produces a covariance matrix
values. The nodal correction and astronomical argu- and correlation coefficients (see Cherniawsky et al.
ment parameters f, u, and V are also embedded in the 2001) that allow a direct method for evaluating the
least squares matrix. (Differences between the old and (in)dependence of the chosen constituents. It is, there-
new approaches are illustrated schematically in Fig. 1.) fore, relatively easy to perform a series of tests to de-
This not only eliminates the need for postfit corrections, termine the best choice of constituents. This issue will
but it also removes the restriction that analysis periods be discussed in more detail later.
should not be much longer than one year (PBL02). That The mathematics underlying multiple inferences is
said, times series longer than 18.6 yr can avoid nodal relatively straightforward. Assume that Ak0 , gk0 , fk0(tj),
corrections completely and are better analyzed using uk0(tj), and Vk0(tj) are the amplitude, phase lag, nodal
techniques that include the satellite constituents directly amplitude correction, nodal phase correction, and as-
(Foreman and Neufeld 1991). tronomical argument, respectively, for the reference
Setting constituent at time tj , whereas Ai , gi , fi(tj), ui(tj), and
Vi(tj) are the analogous values for the i 5 1,Nk0 con-
X k 5 Ak cos gk , stituents to be inferred from that reference constituent.
(5)
Y k 5 Ak sin gk , The nodal correction values and astronomical argu-
ments can be calculated for each time tj , whereas the
Eq. (4) can be re-expressed as the system of j 5 1,m amplitudes and phases can be computed from the har-
linear equations monic analysis. In particular, once the amplitude ratios
n
ri 5 Ai /Ak0 and phase differences ui 5 gk0 2 gi between
the reference and inferred constituents are specified
h(tj ) 5 Z0 1 atj 1 å f k(tj) fX kcos[V k(tj) 1 uk (tj )]
k51 (usually from previous analyses at the same or nearby
1 Y k sin [V k (tj ) 1 uk (tj )]g 1 R(tj ). (6) locations), the only remaining unknowns are Ak0 and
gk0, which can be determined as follows. Setting
If m . 2(n 1 1), this system is overdetermined and is Ck0 j 5 cos [V k0 (tj ) 1 uk0 (tj )],
usually solved by minimizing
Sk0 j 5 sin[V k0 (tj ) 1 uk0 (tj )],
m (8)
Cij 5 cos[V i (tj ) 1 ui (tj )],
å
j51
2
R (tj ) (7) Sij 5 sin[V i (tj ) 1 ui (tj )],

the contribution from the reference constituent and

with respect to the unknowns Z0 , a, and Xk , Yk , for k 5
those to be inferred at time tj,
1,n. Though there are a variety of techniques for per-
forming this least squares fit (F77 used the Cholesky f k0 (tj )Ak0 cos [V k0 (tj ) 1 uk0 (tj )  gk0 ]
algorithm), the approach chosen here is the SVD algo- N
(9)
rithm (Golub and Van Loan 1983; Press et al. 1992) 1 å f i (tj )Ai cos [V i (tj ) 1 ui (tj )  gi ]
described in Cherniawsky et al. (2001). In addition to i51
810 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

can be re-expressed [in analogy with (6)] as

2 3
N
4 5
X k0 f k0 (tj )Ck0 j 1 å
i51
f i (tj )ri (Ci j cos fi  Si j sin fi )

2 3
N
4 5
1 Y k0 f k0 (tj )Sk0 j 1 å
i51
f i (tj )ri (Ci j sin fi 1 Si j cos fi ) ,

(10)
where the unknowns to be determined by solving an
overdetermined system of linear equations are
X k0 5 Ak0 cos gk0 ,
(11)
Y k0 5 Ak0 sin gk0 .
The amplitude and phase of the reference constituent
are then recovered as
qffiffiffiffiffiffiffiffiffiffi
Ak0 5 (X 2k0 1 Y 2k0 ),
(12) FIG. 2. Std CTD lines off southwest Vancouver Island and
gk0 5 tan1 (Y k0 /X k0 ) , western Washington State. Tofino, the tide gauge whose har-
monics were used for all the synthetic tests, and Port Renfrew, the
while those for the inferred constituents are computed reference station that was used to provide inference parameters for
the CTD analysis, are also shown. LA and LB denote sampling
using the prescribed amplitude ratios and phase differ-
lines A and B, respectively.
ences. By simply replacing the coefficients of Xk and Yk
in (6) with those for Xk0 and Yk0 in (10), the inferences
can be included directly in the least squares calculation noise level—as described above. (Yearly analyses of
rather than as postfit corrections. Tofino observations typically have noise-to-signal ratios
of about 17%.) As the period required to separate the
two constituent pairs K1 and P1 and S2 and K2 is ap-
3. Testing the improvements
proximately six months, |v2 2 v1|T ’ 0.33 and the
To assess accuracy, the new methodology was tested Munk and Hasselmann (1964) criterion suggests that
with several synthetic time series that were generated ‘‘meaningful’’ results should be possible by analyzing
using specified constituent amplitudes and phases and both time series without any inference. The same 51
random background noise levels. Constituent nodal constituents selected with the Rayleigh criterion value
corrections and astronomical arguments were computed of 0.33 in F77 were included in all analyses. Table 1 shows
and incorporated at each time level in these syntheses, essentially no difference between the analysis results
consistent with the new analysis approach and what arising from F77 and the new methodology. This indi-
would arise with actual observations. The noise was cates that for relatively short periods such as two months,
created by scaling uniform random numbers that were the accuracy improvements arising from embedding the
generated in the range [20.5, 0.5] with subroutine nodal corrections into the least squares matrix and
RAN1 (Press et al. 1992), so that the ratio of their evaluating the astronomical argument exactly (as op-
standard deviation to the standard deviation of the tidal posed to using a linear approximation) are not appre-
signal was at a prescribed level. ciable. With no background noise, both techniques re-
For the first test, hourly elevations for Tofino, British covered the K1, P1, S2, and K2 signals to within machine
Columbia (Fig. 2), were generated between 1 January precision; however, with 25% noise, they both had ap-
2008 and 1 March 2008 using the same mean sea level proximately 20% errors in the P1 and K2 amplitudes.
and the same amplitudes and phases for the largest eight The matrix condition number for the new methodology,
constituents (Q1, O1, P1, K1, N2, M2, S2, and K2) as are computed as the ratio of the largest to the smallest
employed by the Canadian Hydrographic Service in singular values, was 15.3. So, in this case, the right-hand
their annual tide table predictions (available online side of (3) was a generous upper bound.
at http://www-sci.pac.dfo-mpo.gc.ca/tides_e.htm). Two As a second test, we repeated the previous analysis
time series were computed—one without any back- for only the January portion of the synthesized record
ground noise and the other with a 25% background and now inferred P1 and K2 using their exact amplitude
APRIL 2009 FOREMAN ET AL. 811

TABLE 1. Constituents K1, P1, S2, and K2 true amplitudes and phase lags and their errors when analyzing the 1 Jan–1 Mar 2008 Tofino
synthesized hourly time series with the old (F77) and new harmonic analysis programs. There was no inference, and the noise was uniform
random, with the percentage referring to its std dev relative to the tidal signal std dev.

True Orig—true New—true Orig—true New—true Orig—true New—true Orig—true New—true

amp amp amp amp amp True phase phase phase phase phase
Constituents (m) (no noise) (no noise) (25% noise) (25% noise) (8, UTC) (no noise) (no noise) (25% noise) (25% noise)
P1 0.123 0.000 0.000 0.027 0.026 119.3 20.4 0.0 23.3 23.1
K1 0.388 0.000 0.000 0.020 0.019 121.7 20.1 0.0 22.8 22.2
S2 0.280 0.000 0.000 20.005 20.006 31.5 0.1 0.0 24.5 24.6
K2 0.077 20.001 0.000 0.014 0.014 23.2 0.1 0.0 21.2 21.4

ratio and phase difference relationships relative to K1 leakage in the no-noise case if the inference parameters
and S2, respectively. Table 2 shows analysis results for were not exact, but results from the 25% noise case
the eight constituents included in the synthesis as well as suggest that it would be less than that for F77.
those for NO1, J1, L2, and ETA2, lesser constituents To demonstrate the advantage of embedding the
whose frequencies are close to those inferred. Again, nodal corrections and astronomical arguments in the
the overall accuracy with which both methods recov- least squares matrix, we analyzed in a third test two 12-
ered the amplitudes and phases of the original eight yr time series for Tofino generated with the same eight
constituents is generally quite good. For example, the P1 major constituents and the same noise levels. Though
amplitudes are within 4% of their true values in the assumptions in the nodal correction technique devel-
25% noise case. However, a notable deficiency of the oped by Godin (1972) and employed in F77 become
old inference method is apparent in the amplitudes of increasing invalid as the time series length extends be-
NO1, J1, L2, and ETA2 for the no-noise case. The fact yond 1 yr, a 12-yr analysis can be done with the old
that P1 and K2 energies are present in the time series but method if appropriate array sizes are increased. The
not accounted for in the first stage (least squares fit) of amplitudes and phases arising from these analyses are
the F77 harmonic analysis, causes a leakage to neigh- given in Table 3. Though the phase errors for all eight
boring constituents that is not corrected in the subse- constituents are seen to be reasonably close to zero for
quent postfit inference calculation. In fact, this leakage both methods and for both noise levels, the F77 am-
is still evident in the 25% noise results, as the F77 am- plitude errors for constituents O1, K1, M2, and K2 are
plitudes for these four lesser constituents are consis- approximately 8%, 5%, 2%, and 15%, respectively,
tently larger than those for the new method (whose whereas the analogous values for the new method are
errors arise solely from the background noise). Granted essentially zero, even with the 25% noise-level case.
the new inference method would also display some This improvement is solely a result of having more

TABLE 2. True amplitudes and phase lags and their errors when analyzing the January 2008 Tofino synthesized hourly time series with
the old (F77) and new harmonic analysis programs. Constituents P1 and K2 were inferred from K1 and S2, respectively, using exact
amplitude ratios and phase differences. The noise was uniform random, with the percentage referring to its std dev relative to the tidal
signal std dev.

True Orig—true New—true Orig—true New—true True Orig—true New—true Orig—true New—true
amp amp amp amp amp phase phase phase phase phase
Constituents (m) (no noise) (no noise) (25% noise) (25% noise) (8, UTC) (no noise) (no noise) (25% noise) (25% noise)
Q1 0.044 20.001 0.000 20.001 20.001 113.8 0.3 0.0 4.4 4.2
O1 0.247 0.004 0.000 0.004 0.001 116.5 20.5 0.0 0.9 1.3
NO1 0.0 0.019 0.000 0.019 0.016 0.0 — — — —
P1 0.123 0.000 0.000 20.005 20.005 119.3 0.1 0.0 2.8 2.6
K1 0.388 0.001 0.000 20.016 20.017 121.7 0.1 0.0 2.8 2.6
J1 0.0 0.011 0.000 0.014 0.007 0.0 — — — —
N2 0.204 20.001 0.000 0.008 0.010 349.8 0.3 0.0 3.3 3.0
M2 0.986 0.002 0.000 0.020 0.017 10.0 20.3 0.0 20.3 0.0
L2 0.0 0.014 0.000 0.019 0.015 0.0 — — — —
S2 0.280 0.0 0.000 0.008 0.007 31.5 0.6 0.0 0.8 0.2
K2 0.077 0.0 0.000 0.002 0.002 23.2 0.6 0.0 0.8 0.2
ETA2 0.0 0.006 0.000 0.013 0.007 0.0 — — — —
812 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

TABLE 3. As in Table 2 except when analyzing the 1 Jan 2000–31 Dec 2011 Tofino synthesized hourly time series with the old (F77) and
new harmonic analysis programs. The noise was uniform random, with the percentage referring to its std dev relative to the tidal signal sd
dev.

True Orig—true New—true Orig—true New—true True Orig—true New—true Orig—true New—true
amp amp amp amp amp phase phase phase phase phase
Constituents (m) (no noise) (no noise) (25% noise) (25% noise) (8, UTC) (no noise) (no noise) (25% noise) (25% noise)
Q1 0.044 20.003 0.000 20.003 20.001 113.8 2.4 0.0 1.6 20.8
O1 0.247 20.020 0.000 20.021 20.001 116.5 1.0 0.0 0.7 20.3
P1 0.123 0.001 0.000 20.001 20.005 119.3 0.0 0.0 0.3 0.3
K1 0.388 20.021 0.000 20.022 20.001 121.7 20.6 0.0 20.4 0.2
N2 0.204 20.004 0.000 0.004 0.000 349.8 20.5 0.0 20.4 0.0
M2 0.986 0.021 0.000 0.022 20.001 10.0 20.4 0.0 20.4 0.0
S2 0.280 20.001 0.000 0.000 0.001 31.5 0.0 0.0 20.1 20.1
K2 0.077 20.012 0.000 20.012 20.001 23.2 21.8 0.0 22.9 20.8

accurate nodal corrections in the new method. The method rather than the Foreman and Neufeld (1991)
reason the phases were relatively consistent between method for time series fewer than 18.6 yr.
methods is that the analysis period was approximately
centered over 2006, which corresponds to K1/O1 maxi-
4. Constituent selection and error estimates
mum amplitudes (and M2/N2 minimum amplitudes) in
their 18.6-yr variations, thereby making u(t0) [as defined To illustrate how SVD-generated correlations and
in Eq. (1)] a good approximation. However, the same error estimates can be used to assist in constituent se-
centering means that the associated amplitude correc- lection, we first return to our analysis of the January
tions f(t0) for the old method were too large for the 2008 portion of the synthetic Tofino time series with a
diurnals and too small for M2 and N2. 25% noise-to-signal ratio. This analysis used the same
The foregoing 12-yr analysis was also carried out with 38 constituents that were chosen with a Rayleigh con-
the Foreman and Neufeld (1991) long (18.6 yr and stant of 0.97 and the Rayleigh criterion decision tree
more) analysis software, which avoids nodal corrections described in F77. Constituents P1 and K2 were inferred
completely by directly including both satellites and from K1 and S2, respectively, using the exact amplitude
major constituents in the least squares fit. A standard ratios and phase differences used in the synthesis. The
run with the full complement of satellite and major RMS deviation for the original time series was 0.838 m,
constituents solves for 529 pairs of amplitudes and pha- while the RMS residual after the least squares fit was
ses, and because many of these have nearest neighbors 0.200 m, a value consistent with our expected 25% noise
with a frequency separation of 18.621 year 21, the matrix level.
condition number can be large when the analyzed time As described in Cherniawsky et al. (2001), the SVD
series length is much less than 18.6 yr. —such is the case least squares solution also produces estimates of the
here. With essentially no background noise (it only arises correlation between all possible Xk and Y [Eq. (5)]
from rounding the synthesized hourly values to three constituent combinations and error estimates for these
decimal places), there were errors of 4%, 5%, and 13% same variables. However, these error estimates assume
in the O1, K1, and M2 amplitudes, respectively (Table 3). a normal distribution (Press et al. 1992) for the R(tj)
Increasing the resolution to four decimal places scarcely residuals of Eq. (6), and as pointed out by Munk et al.
changed the results, so the matrix condition number must (1965), PBL02, and others, the residual spectrum is
be the term dictating accuracy in Eq. (3). Reducing the generally more red than white, with cusps around each
number of constituents/satellites in the analysis to only tidal frequency. So, the assumption of a normal distri-
those included in the synthesis improved the accuracy bution is questionable. If the time series had a constant
slightly (presumably by reducing the matrix condition sampling interval, then it would be relatively easy to
number), while increasing the time series length to 19 yr follow PBL02 and use fast Fourier transform (FFT)
produced the same accuracy as that for the new method, methods to estimate the background variance in a fre-
even with the full complement of 529 constituents. These quency band around each—or at least each ma-
tests demonstrate that if there is no reason to question jor—constituent, and then apply a parametric bootstrap
the satellite amplitude ratios and phase differences based method to provide better uncertainty estimates. How-
on tidal potential theory and inherent in standard nodal ever, with the provision for irregular sampling within
corrections, it is more accurate to use the new analysis the analyzed time series, it is not obvious how a single
APRIL 2009 FOREMAN ET AL. 813

approach can be used to compute background variance

constituents—
values with 12
TABLE 4. As in Table 2 except when analyzing randomly selected values within the periods of 1 Jan–30 Jun 2008 (columns 3, 4, 8, 9) and 1 Jan–1 Mar 2008 (columns 5, 6, 10, 11).

true phase
Phase 488
estimates around tidal frequencies. Standard FFT

1.7
21.0

1.5
1.5

0.9
20.8

1.2
1.3
—

—
techniques can be used with regular sampling but they
cannot be used with irregular. Perhaps a more general
Fourier approach that employs a least squares approach
to find the energy at specific frequencies separated by
constituents—
values with 38

true phase
Phase 488

22.2 regular intervals could be developed, but such an ap-

20.3

4.3
2.2

0.5
20.7

7.7
22.8
proach raises issues of what the interval and frequencies
—

—
should be. Alternatively, if a long regularly sampled time
series is available at the same or a nearby location to the
one being analyzed, then a de-tided spectra and more
values with 11
constituents—

accurate amplitude and phase estimates could be com-

true phase
Phase 744

10.7
22.8
1.1

20.4
20.3

21.2
0.0

1.0
puted following a bootstrap approach like that described
—

—
in PBL02. Although this is an intriguing problem, it is
beyond the scope of the present study. So for now, we are
left with the estimates presently provided by an admit-
values with 51
constituents—

tedly incorrect Gaussian assumption for the residuals.

true phase
Phase 744

For the January 2008 analysis of synthetic Tofino el-

11.6
24.8
1.2

22.4
0.2

22.0
20.2

1.7
—

evations, the largest correlations were 0.157 between

the Yk component of ETA2 and the Xk component of S2,
the ratio of the largest to the smallest singular value was
2.11, while the error estimates for all the Xk and Yi
(8, UTC)
phase

10.0

31.5
23.2
0.0

0.0

0.0

0.0
113.8
116.5

119.3
121.7

349.8
True

constituent amplitudes ranged between 0.006 and 0.011

m. Despite the incorrect Gaussian assumption under-
lying their calculation, the actual amplitude errors listed
values with 12

in the sixth column of Table 2 compare favorably to

constituents
—true amp
Amp: 488

0.015
0.008

20.004
20.012
0.041
0.007
0.025

0.008
0.002

these Xk and Yi error estimates. Apart from Z0, the

—

inferred constituents P1 and K2, and the remaining six

constituents used in the synthesis (Q1, O1, K1, N2, M2,
and S2), only MU2, MK3, SK3, S4, 2SM6, and M8 had
amplitudes that were more that twice as large as their
constituents—
values with 38
Amp: 488

true amp

standard deviation estimates. Thus, if we could assume

0.022
0.008
0.020
20.012
20.007
0.038
0.006
0.019
0.009
0.021
20.030
0.004

that the postfit residuals [R(tj) in Eq. (1)] were Gaus-

sian, a Student’s t test would only find these constituents
to be statistically different from zero at the 95% level.
Repeating the analysis with only these ‘‘significant’’
constitituents—
values with 11

constituents plus Z0, and inferring P1 and K2, reduced

Amp: 744

true amp
20.002
0.004
0.047
0.023
20.013

20.012
0.002
0.020
20.010
20.013

the RMS postfit residual slightly to 0.198 m. It also de-

—

creased the largest correlations (now between M2 and

N2) to 0.105, reduced the range of the Xk and Yk error
estimates to 0.008 from 0.11 m, and maintained essen-
tially the same accuracy (amplitudes to within 1 mm and
values with 51
constituents
—true amp
Amp: 744

phases to within 0.18) as that shown in the sixth column

20.002
0.007
0.045
0.022
20.009
0.005
20.011
0.003
0.023
20.014
20.011
0.010

of Table 2; thereby suggesting that the additional 23

constituents did not really contribute to the initial
analysis.
The previous test used hourly sampled data for which
0.044
0.247

0.123
0.388

0.204
0.986

0.280
0.077
True
amp
(m)

0.0

0.0

0.0

0.0

the Rayleigh (F77), Munk and Hasselmann (1964), or

Foreman and Henry (1989) criteria could provide rea-
Constituents

sonable guidance on constituent selection. In this next

series of tests, we analyze randomly sampled data, so
ETA2

that the first two criteria do not apply. The dataset is

NO1

M2
Q1
O1

K1

N2

K2
L2
P1

S2

based on the same 2008 synthesized Tofino time series

J1
814 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

used previously but in this case, we randomly select of Vancouver Island (Fig. 2). Eleven and 12 respective
which hourly samples are to be used in the analysis. In stations along these lines have generally been sampled
the first test, we used the Press et al. (1992) subroutine 2–3 times a year since 1980, with observations taken at
RAN1 to randomly choose 744 hourly ‘‘observations’’ the standard depths of 0, 5, 10, and every 10 m thereafter
(with 25% noise) during the 1 January–30 June period. down to the bottom, or 2400 m. Though the stratifica-
(As an aside, these observations need not be ordered tion and internal tide patterns do change seasonally, it is
chronologically for the analysis.) Table 4 shows the feasible to restrict each CTD time series to one season
amplitudes and phases obtained by analyzing with the and analyze for tidal variations in salinity and temper-
same set of 51 constituents that would be selected by the ature. Here we restrict the observations to June through
F77 decision tree for a 6-month record. Not only are the September, inclusive. There are between 18 and 52
results reasonably accurate but the ratio of maximum to observations at each standard depth at each station, and
minimum singular value (matrix condition number) was the analyses solve directly for a linear trend and con-
2.78 and the maximum correlation coefficient was 0.137, stituents Z0, K1, M2, while inferring Q1, O1, and P1 from
between the lesser constituents OO1 and UPS1. Re- K1, and N2, S2, and K2 from M2. Inference parameters
moving all those constituents with amplitudes less than were taken from a 1-yr analysis of hourly tide gauge
twice the error estimates and rerunning the analysis observations at Port Renfrew (Fig. 2). The objectives of
reduced the analysis set to only 20 constituents. The the analyses are to (i) compute the magnitude of the
matrix condition number was now 2.40 and the largest tidal variations in these observations; (ii) determine if
correlation coefficient was 0.115, between UPS1 and the M2 variations show evidence of internal tides (or at
TAU1. As shown in Table 4, the amplitude accuracy for least that portion of the internal tide that is phase locked
this new analysis is close to the previous one with 51 with the barotropic tide); and (iii) determine if the ob-
constituents and with one exception (K1): the phase servations show a linear trend.
errors have decreased. The linear trend in temperature and the M2 tempera-
The next tests analyzed 488 randomly selected hourly ture phases along line B are shown in Figs. 3 and 4, re-
values during the period of 1 January–1 March with the spectively. Average summer seasonal currents crossing
same constituent dataset that F77 would choose for that this line include a near-surface shelf break current (SBC)
period, plus P1 and K2. The matrix condition number flowing to the southeast, a near-surface Vancouver Island
was now 4.94 and as would be expected, the largest coastal current (VICC) flowing to the northwest, and a
correlation coefficients ranged between 0.822 and 0.832 California Undercurrent (CUC) flowing to the northwest
for P1/K1 and K2/S2. Amplitude and phase errors for the along the continental slope and centered at about 200-m
major constituents are shown in Table 4. Repeating the depth (Freeland et al. 1984; Foreman et al. 2000). A
analysis with P1 and K2 inferred and eliminating all warming of up to 0.058C yr21 in the SBC, a cooling of up
constituents with amplitudes less than twice their error to 0.038C yr21 in the VICC, and no change in the CUC
estimates improved the accuracy (Table 4) of all con- are evident in Fig. 3. However, the standard deviation
stituents involved in the inference, with the exception of estimates associated with these values generally range
the K1 amplitudes. The matrix condition number between 0.018 and 0.038C yr21, so these trends are not
dropped to 1.71, and the largest correlation coefficient statistically significant. Correlation coefficients between
was now 0.155, between Q1 and O1. a and the other constituent parameters are generally
Though many more tests could be performed, the less than 0.3.
preceding few have demonstrated that by monitoring Though the M2 phase patterns seen in Fig. 4 are noisy,
correlation coefficients, matrix condition numbers, and there is the suggestion of a vertical mode structure at the
Student’s t test values (albeit based on the generally edge of the continental shelf where internal tides are
incorrect assumption that the residuals have a Gaussian known to be generated (Drakopoulos and Marsden
distribution), an iterative procedure can be used to de- 1993). However, the pattern might also be due to bar-
termine the best set of constituents to be included with otropic tide advection of a temperature field with ver-
this new harmonic analysis. tical and lateral structure. As with the linear trend re-
sults, the relatively few number of points in the time
series means that these particular results are not sta-
5. New applications tistically significant. Nevertheless, it has been demon-
strated that long-term CTD analyses are feasible with
a. CTD analyses
this new program so that as the time series continues
The first example is the analysis of CTD observations to lengthen, statistically significant results might be ex-
along lines A (LA) and B (LB) off the southwest coast pected.
APRIL 2009 FOREMAN ET AL. 815

FIG. 4. The M2 temperature phase lags (8) along line B.

FIG. 3. The linear trend in temperature (8C yr21) along line B and
locations of VICC and SBC.

crosses the Strait of Georgia in a southeasterly direction

It should be mentioned that these randomly sampled (Fig. 5). These points are separated by approximately
CTD time series could have been analyzed with the 5.7 km. The fact that the strait is only about 30 km wide
software that is described in Foreman and Henry (1979, means that there is frequent data loss as a result of
henceforth FH79), which employs the same astronom- signal contamination from nearby land, particularly at
ical, nodal, and inference correction approach as F77. the northern and southern portions of the track. Out of
However, for analyses of relatively few observations a maximum sample size of 364, there are only 12 loca-
spanning 25 yr, only a few constituents can be resolved. tions with 58 or more noncontaminated values and the
The FH79 approach should be less accurate than the best site has only 276. Harmonic analyses were per-
new approach, because it assumes constant nodal cor- formed at each of these 12 locations using 17 constitu-
rection parameters over the duration of the analysis pe- ents composed of Z0, Sa, Ssa, and the 7 largest diurnal
riod and it does not allow multiple inferences. Though and 7 largest semidiurnal constituents. No inference was
the true temperature variations at tidal amplitudes are performed initially, though it is well-known (Parke et al.
not known in this case, a simple analysis of the 5-m 1987; Ray 1998; Cherniawsky et al. 2001) that the
temperature data at the line B station nearest to the 9.9156-day sampling interval for T/P can cause signifi-
shore seems to confirm this supposition. Using FH79 to cant aliasing, even with records as long as 10 yr. Note-
solve directly for constituents O1, K1, M2, and S2 and worthy constituent pairs that may remain difficult to
inferring Q1, P1, N2, and S2, respectively, from those four separate with this record length are K1 and Ssa and P1
produced M2 and K1 amplitudes that were 20% and 17% and K2, and this was borne out in the correlation coef-
larger, respectively, and phases that were 178 smaller and ficients arising from the analyses. For the central loca-
518 larger, respectively, than those described above with tion with the sample size of 276, the largest correlation
the new method. The fact that the overall K1 differences coefficients were 0.270 and 0.226 between the Xk and Yk
(FH79 minus new method) are larger than those for M2, components of K2 and P1, respectively, while the next
while the M2 temperature amplitudes themselves are largest values were 0.193 and 0.146 for the analogous
larger, is consistent with what would be expected to be components of K1 and Ssa, respectively. At a site with a
larger errors in the FH79 nodal corrections during that sample size of only 119 values, the K2/P1 values rose to
period. 0.340 and 0.290. Figure 6 illustrates the relationship
between the along-track P1 and K2 amplitudes, their
b. Satellite altimeter analyses
correlation coefficients, and the analysis sample sizes.
The second example is the analysis of altimeter data (The correlation coefficients at sites 25 and 36 are not
from the TOPEX/Poseidon satellite. In this case, the shown, because these analyses had numerous constitu-
time series are restricted to the period from 30 Sep- ent separation problems, thereby making all results
tember 1992 to 8 August 2002 and are taken from col- questionable.) Here, P1 and K2 amplitudes at the
located points along a small portion of track 90 that nearby Point Atkinson tide gauge (Fig. 5) are also
816 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 26

FIG. 6. Select harmonic analysis results for P1 and K2 at collocated

sites along the T/P ground track in the Strait of Georgia.
FIG. 5. T/P satellite ground track 90 and collocated stations in
the Strait of Georgia. The location of the reference tide gauge at
Point Atkinson is also shown.
with and without background noise. Where feasible,
comparisons were also done with results from F77 and
shown as a means of evaluating accuracy. Clearly, the the Foreman and Neufeld (1991) long analysis codes.
smaller sample sizes lead to higher correlations and The use of correlation matrix output in constituent se-
generally less accurate results. Repeating these analyses lection decisions was demonstrated with two further
with K2 inferred from S2, and the inference relationships examples. Finally, the software was applied to two
computed from the Point Atkinson harmonics, im- problems that could not have been solved with the older
proved the P1 and K2 accuracy at the sites near the ends software. The first application was the analysis of 25 yr
of the track. For example, the K2 amplitude at site 36 of CTD data along a transect off southwest Vancouver
decreased from 133.6 to 68.3 cm, much closer to the Island that suggested both long-term temperature
61.9-cm value at Point Atkinson. trends in the shelf break and Vancouver Island Coastal
As with the CTD analysis example, this application Currents as well as vertical variations in M2 phase that
could also have been solved with the FH79 software. might be attributable to internal tides. The second ap-
However, for the same reasons explained earlier, we plication was the analysis of T/P satellite altimetry data
would expect poorer accuracy, because the nodal cor- along a ground track in the Strait of Georgia whose
rection parameters would be held constant during the proximity to land has led to significant data dropout. In
10-yr period of the analysis. this case, the relationships between the number of valid
data points and constituent aliasing is clearly seen in the
correlation matrices and thus can be used as a guide to
6. Summary and conclusions
determine which constituents should be included di-
The previous presentation has described and applied rectly in the analysis and which should be inferred.
new computer software that permits more versatility in Though the set of astronomical constants and con-
the harmonic analysis of tidal time series. Specific im- stituents used in this new software is the same as that
provements to traditional methods include the analysis employed in F77 and PBL02—namely, those derived
of randomly sampled and/or multiyear data; more ac- from Cartwright and Tayler (1971) and Cartwright and
curate nodal correction, inference, and astronomical Edden (1973)—they could easily be replaced with more
argument adjustments through direct incorporation into recent versions such as those of Hartmann and Wenzel
the least squares matrix; and correlation matrices and (1995; available online at http://bowie.gsfc.nasa.gov/
error estimates that facilitate decisions on the selection hw95/), which also include the tide generating potential
of constituents for the analysis. One- and two-dimensional of the planets Venus, Jupiter, Mars, Mercury, and Sat-
time series can be analyzed with the same code and in the urn. Only changes to the astronomical input file and list
case of the latter, the final harmonics are also expressed in of constituents to be included in the analysis would be
terms of current ellipse parameters. required.
The accuracy of the new methodology was assessed This new software (in FORTRAN) is freely available as
through a series of test analyses using synthetic data part of the Institute of Ocean Sciences (IOS) Tidal
 APRIL 2009 FOREMAN ET AL. 817

Package and can be downloaded, along with sample input ——, and E. M. Neufeld, 1991: Harmonic tidal analyses of long time
data and a short explanatory readme file, from the series. Int. Hydrogr. Rev., LXVIII, 85–108.
——, R. E. Thomson, and C. L. Smith, 2000: Seasonal current
FTP site given on http://www-sci.pac.dfo-mpo.gc.ca/osap/
simulations for the western continental margin of Vancouver
projects/tidpack/tidpack_e.htm. All comments and prob- Island. J. Geophys. Res., 105, 19 665–19 698.
lem reports are welcome. Freeland, H. J., W. R. Crawford, and R. E. Thomson, 1984: Cur-
rents along the Pacific coast of Canada. Atmos.–Ocean, 22,
Acknowledgments. We thank Trish Kimber, Wendy 151–172.
Callendar, and Ming Guo for their assistance with the Godin, G., 1972: The Analysis of Tides. University of Toronto
figures; Jake Galbraith for processing the CTD data; Brian Press, 264 pp.
Golub, G. H., and C. F. Van Loan, 1983: Matrix Computations.
Beckley and Richard Ray for providing the T/P altimeter
The Johns Hopkins University Press, 476 pp.
data; and two anonymous reviewers for their constructive Hartmann, T., and H.-G. Wenzel, 1995: The HW95 tidal potential
comments on an earlier version of the manuscript. catalogue. Geophys. Res. Lett., 22, 3553–3556.
Jay, D. A., and E. P. Flinchem, 1997: Interaction of fluctuating
REFERENCES river flow with a barotropic tide: A test of wavelet tidal
analysis methods. J. Geophys. Res., 102, 5705–5720.
Cartwright, D. E., 1999: Tides: A Scientific History. Cambridge ——, and ——, 1999: A comparison of methods for analysis of
University Press, 292 pp. tidal records with multi-scale non-tidal background energy.
——, and R. J. Tayler, 1971: New computations of the tide-gen- Cont. Shelf Res., 19, 1695–1732.
erating potential. Geophys. J. Roy. Astron. Soc., 23, 45–74. Munk, W., and K. Hasselmann, 1964: Super-resolution of tides.
——, and C. A. Edden, 1973: Corrected tables of tidal harmonics. Studies on Oceanography—A Collection of Papers Dedicated to
Geophys. J. Roy. Astron. Soc., 33, 253–264. Koji Hidaka, K. Yoshida, Ed., University of Tokyo, 339–344.
Cherniawsky, J. Y., M. G. G. Foreman, W. R. Crawford, and R. F. ——, B. Zetler, and G. W. Groves, 1965: Tidal cusps. Geophys. J.
Henry, 2001: Ocean tides from TOPEX/Poseidon sea level Int., 10, 211–219.
data. J. Atmos. Oceanic Technol., 18, 649–664. Ortega, J. M., 1972: Numerical Analysis: A Second Course. Aca-
Doodson, A. T., 1921: The harmonic development of the tide demic Press, 201 pp.
generating potential. Proc. Roy. Soc. London, A100, 306–328. Parke, M. E., R. H. Stewart, D. L. Farliss, and D. R. Cartwright,
Drakopoulos, P. G., and R. F. Marsden, 1993: The internal tide off 1987: On the choice of orbits for an altimetric satellite to study
the west coast of Vancouver Island. J. Phys. Oceanogr., 23, ocean circulation and tides. J. Geophys. Res., 92, 11 693–11 707.
758–775. Parker, B. B., 2007: Tidal analysis and prediction. NOAA Special
Foreman, M. G. G., 1977: Manual for tidal heights analysis and Publication NOS CO-OPS 3, U.S. Department of Commerce,
prediction. Pacific Marine Science Rep. 77-10, Institute of 378 pp. [Available online at http://tidesandcurrents.noaa.gov.]
Ocean Sciences, Patricia Bay, 66 pp. [Available online at http:// Pawlowicz, R., B. Beardsley, and S. Lentz, 2002: Classical tidal
www.pac.dfo-mpo.gc.ca/SCI/osap/publ/online/heights.pdf.] harmonic analysis including error estimates in MATLAB
——, and R. F. Henry, 1979: Tidal analysis based on high and low using T_TIDE. Comput. Geosci., 28, 929–937.
water observations. Pacific Marine Science Rep. 79-15, Institute Press, W. H., S. A. Teukolsky, W. T. Vettering, and B. P. Flannery,
of Ocean Sciences, Patricia Bay. [Available online at http:// 1992: Numerical Recipes in FORTRAN. 2nd ed. Cambridge
www.pac.dfo-mpo.gc.ca/SCI/osap/publ/online/high-low.pdf.] University Press, 963 pp.
——, and ——, 1989: The harmonic analysis of tidal model time Ray, R. D., 1998: Spectral analysis of highly aliased sea-level sig-
series. Adv. Water Resour., 12, 109–120. nals. J. Geophys. Res., 103, 24 991–25 003.

View publication stats