Ocean Tides Modeling using

Satellite Altimetry

by

Hok Sum Fok

Report No. 501

Geodetic Science
The Ohio State University
Columbus, Ohio 43210
December 2012
Ocean Tides Modeling using Satellite Altimetry

by
Hok Sum Fok

Report No. 501

Geodetic Science
Ohio State University
Columbus, Ohio 43210

December 2012
 Preface

This report was prepared for and submitted to the Graduate School of the Ohio State
University as a dissertation for partial fulfillment of the requirements for the Doctor of
Philosophy (PhD) degree.

This research is conducted under the supervision of Professor C.K. Shum, Division of
Geodetic Science, School of Earth Sciences, The Ohio State University. This research is
partially supported by grants from NASA’s Physical Oceanography program, including
Ocean Surface Topography Science Team Program (JPL 1356532, JPL 1384376, and
UC154-5322), and from Ohio State University's Climate, Water, and Carbon (CWC)
Program (http://cwc.osu.edu).

i
 Abstract

Ocean tides, resulting from the gravitational attractions of the Moon and the Sun,
represent 80% of the ocean surface topography variability with a practical importance for
commerce and science over hundreds of years. Tides have strong influence on the
modeling of coastal or continental shelf circulations, play a significant role in climate due
to its complex interactions between ocean, atmosphere, and sea ice, dissipate their energy
in the ocean and solid Earth, and decelerate the Moon’s mean motion. Oceanographic
studies and applications, including coastal or continental shelf ocean circulations, also
require observations to be ‘de-tided’ using ocean tidal forward prediction models before
geophysical or oceanographic interpretation, particularly over coastal regions. Advances
in satellite radar altimetry technology enabled a globally sampled record of sea surface
height (SSH) and its changes over the past two decades, particularly after the launch of
TOPEX/POSEIDON satellite. This geophysical record enables numerous scientific
studies or discoveries, including improved global ocean tide modeling.

Several contemporary ocean tide models have been determined either through the
assimilation of satellite altimetry and coastal tide gauge data, often referred to as
‘assimilation models’ (e.g. FES2004, NAO.99b and TPXO6.2/7.1/7.2), or via the use of
altimetry observations in an ‘empirical modeling’ approach to solve for tidal constituents
based on a-priori tide models, including assimilated models (e.g. DTU10,
EOT08a/10a/11a, GOT00.2/4.7). However, ocean tide model accuracy is still much
worse, up to an order of magnitude, in the coastal regions or over partially or
permanently sea-ice or ice-shelf covered polar ocean, than that of models in the deep
ocean.

Here observation-based, empirical ocean tide models with 0.25o×0.25o spatial resolution,
the OSU12 models, has been determined using improved multi-satellite altimetry data
from TOPEX, Jason-1/-2, Envisat, and GFO, and based on a novel approach via spatio-
temporal combination, along with a robust estimation technique. We first demonstrate the
effectiveness of the spatio-temporal combination approach when comparing with various
ocean tide solutions under different data weighting schemes (i.e. equally weighted
solution, the weighted solution based on spatial (co-)variances, and the weighted solution
based on temporal (co-)variances). The generated tide models show substantial
improvement near coastal regions when compared against contemporary ocean tide
models using assessment from independent tidal constants of tide gauges and from
variance reduction studies using altimetry data. The improvement is particularly apparent
in regions with high hydrodynamic variability, yet the model accuracy is still region-
dependent. The model is available at: http://geodeticscience.org/oceantides/OSU12v1.0/

For the first time, the potential seasonality of ocean tides in subarctic regions has been
demonstrated. A statistically significant difference in variance reduction of multi-mission
altimeter SSH anomalies are observed in the subarctic ocean study region during summer
ii
and winter seasons. The variability in the SSH anomalies during winter are 15–30%
larger than those of summer, and we hypothesize that seasonality of tides contributes to
the observed SSH variability. The subsequent seasonal ocean tide solutions estimated
using observations only in the winter and in the summer seasons, reveal detectable
seasonal tidal patterns in the Chukchi Sea near the eastern Siberia region where it is
known to have seasonal presence of sea-ice covers.

iii
 Table of Contents

Preface.................................................................................................................................. i
Abstract ............................................................................................................................... ii
Table of Contents............................................................................................................... iv
Chapter 1: Introduction ....................................................................................................... 1
1.1 Ocean Tides and its Modeling Development ....................................................... 1
1.2 Contemporary Ocean Tide Models Using Satellite Altimetry ............................. 3
1.3 Problem Statement and Research Methodology .................................................. 5
Chapter 2: Theory of Ocean Tides and Tidal Analysis....................................................... 8
2.1 Tide-generating Potential and its Species ............................................................ 8
2.2 Harmonic Expansion of the Tide-generating Potential ...................................... 11
2.3 Tidal Analysis .................................................................................................... 14
2.3.1 Laplace Tidal Equations ................................................................................... 14
2.3.2 Harmonic Analysis ........................................................................................... 16
2.3.3 Response Analysis............................................................................................ 17
2.3.4 Orthotide Formulation and its Relation to Harmonic Constants ...................... 19
2.4 Chapter Summary............................................................................................... 20
Chapter 3: Satellite Altimetry Technique ......................................................................... 22
3.1 Sea Surface Height measurement in satellite altimetry...................................... 22
3.1.1 Measurement principle of satellite altimetry .................................................... 25
3.1.2 Corrections to altimeter radar measurements ................................................... 28
3.2 Tidal Aliasing..................................................................................................... 32
3.3 Chapter Summary............................................................................................... 35
Chapter 4: Ocean Tides Modeling using Satellite Altimetry............................................ 36
4.1 A Review of Ocean Tides Modeling Approaches Using Satellite Altimetry .... 36
4.1.1 Hydrodynamic modeling .................................................................................. 36
4.1.2 Assimilation modeling...................................................................................... 38
4.1.3 Empirical modeling .......................................................................................... 39
4.2 A Novel Spatio-Temporal Combination Approach ........................................... 42
4.2.1 Observation Equation Setup ............................................................................. 43

iv
 4.2.2 Weighted Least-Squares Solution, Representation of Spatio-Temporal
Combination, and Covariance Model Specification.................................................. 44
4.2.3 Techniques for Spatio-Temporal Combination ................................................ 52
4.2.4 Outlier Detection and Robust Estimation......................................................... 55
4.2.5 Altimetry Data Used and Computational Procedures....................................... 59
4.3 Chapter Summary............................................................................................... 62
Chapter 5: Ocean Tide Modeling Solutions and Accuracy Analysis ............................... 63
5.1 External Accuracy Assessment Using Tide Gauges and Tidal Solutions under
Different Weighting Schemes ....................................................................................... 63
5.2 Internal Accuracy Assessment Using Variance Reduction Test........................ 69
5.3 Global Result and Coastal Assessment .............................................................. 77
5.4 Seasonal Variations of Ocean Tides in the Sub-Arctic Ocean........................... 87
5.5 Chapter Summary............................................................................................... 96
Chapter 6: Conclusion and Future Work .......................................................................... 97
Bibliography ..................................................................................................................... 99
Appendix A: Tables for Overall Results......................................................................... 112

v
 Chapter 1: Introduction

1.1 Ocean Tides and its Modeling Development

Ocean tides, which is a kind of the rise and fall of the sea level as observed from
coasts around the world, are a geophysical phenomenon resulting from the gravitational
attraction of the Sun and the Moon acting on the Earth. They are the principal variable
component of ocean surface topography of the global ocean (Ray, 1993; Le Provost et al.,
1994). They have a strong influence on modeling coastal or continental shelf circulations,
and on short-term (e.g., half daily, daily, fortnightly, monthly) mass variations, with the
later having a significant impact on the Earth rotation and on the change the length of day
(Gross, 1993; Ray et al., 1994). They dissipate their energy in the ocean and solid Earth,
decelerate the Moon’s mean motion, and play a significant role in climate due to its
complex interactions between ocean, atmosphere, and sea ice (Shum et al, 2001).
Historically, one of the practical applications of the ocean tide prediction is global trading
and commerce. Ocean tides have been used for navigation safety when ships approach
the harbors and the coasts.
Due to the fundamental importance of ocean tides, scholars engaged themselves to
advance scientific understanding of ocean tides during the past two centuries. Since the
establishment of Newton’s equilibrium theory which explained the forces that generate
tides, Laplace (1775) formulated his dynamic equations, called Laplace Tidal Equations
(LTE), based on the concept of a hydrodynamic response of the ocean to the tide-
generating force in the form of partial differential equations. Despite the plausible theory,
it is impossible to obtain an analytical solution since the equations strongly depend on the
bathymetry and the shape of the ocean basin. The account for the non-linearity, such as
bottom friction and advection terms, also renders the partial differential equations to be
solved numerically nowadays. Until the late 19th and early 20th century, Darwin (1883)
provided a practical and efficient method, called harmonic analysis, for the tidal analysis 1
and prediction. Newcomb’s solar theory (1895) and Brown’s lunar theory (1905) on the
solar and lunar ephemerides, respectively, facilitated the initial development of tidal
sciences. Doodson (1921) formulated a set of numbers, called Doodson Numbers, to
describe the major tidal frequencies. The aforementioned scholars provided a solid
foundation to ocean tide theory and its practical prediction methodology along with
future development.
With the continuous interest in the physics of the global ocean and its circulation, and
their effect on climate and the Earth system as a whole, sea level was measured by tide
gauges along continental coastlines at 739 coastal sites (refers to coastal tide gauges), and

1
Tidal analysis refers to the usage of the sea level measurements to solve for the amplitude and phase lag
of different tidal constituents, according to their respective astronomical frequencies, via least-squares.
1
by tide gauges or bottom pressure recorders at islands at 102 deep-sea sites (refers to
pelagic tide gauges) 2 . Continuous installation of additional tide gauges is anticipated.
Though the understanding of ocean tides has been further improved based on the long-
term observations, our knowledge of ocean tides remains limited in the vicinity of
coastlines and mid-ocean islands. Moreover, sea level measurements far away from
coasts made by bottom pressure gauges is slow and not cost-effective (Pugh, 1987).
Satellite altimetry, an active remote sensing technique introduced in the 1970s,
provided a synoptic means of observing the global ocean and other surfaces on the Earth
as well as extra terrestrial bodies. On the Earth, this initiative was driven by the initial
quest to improve our understanding of global ocean circulations. Due to its orbital design,
the satellite revisits at the same location for each cycle at a pre-defined temporal period
along its ground tracks. Therefore, the acquired sea surface height (SSH) data from
satellite altimetry can be analogously served as tide gauge measurements at each location
in a global sense. Despite its spatial and global coverage, the temporal sampling is
relatively sparse when compared to tide gauge records.
Thanks to precise orbit determination techniques, along with the improvement in
instrumental, media, and geophysical corrections, the accuracy of satellite altimetry
improved significantly since the launch of TOPEX/POSEIDON in 1992. Owing to its
near-global coverage, the altimeter range accuracy of ± (2–3) cm, and the optimal orbital
configuration for adequate spatial (e.g., equatorial cross-track separation of ~100 km) and
temporal (i.e., approximately weekly) sampling, it makes a significant advancement in
the study of the general circulation of the world’s ocean (Fu et al., 1994). Since then,
other missions, such as ERS-2, GFO, Envisat, Jason-1, and Jason-2, have been launched
to extend the geophysical and oceanographic observations, in particular for the study of
ocean’s role in climate change, including sea-level rise, general ocean circulation and
heat transport.
Before the launch of TOPEX/POSEIDON, the first altimetry-data derived ocean tide
model by Cartwright and Ray (1990a, b, 1991) is based on 2.5 years of altimetry data
from Geosat. It is more accurate than the, then state-of-the-art, hydrodynamic tide model
derived by Schwiderski (1980), which was derived using bathymetry data and
assimilating the tidal constants (i.e., amplitude and phase of tides) computed from global
tide gauge records into a hydrodynamic tidal equation via a finite difference method.
Since then, more than ten ocean tide models were developed in the mid-1990s,
primarily because of the availability of TOPEX/Poseidon altimetry and the fact that
ocean tides are the major contributor to SSH variation. These models serve the
oceanographic and geophysical communities who are ultimately interested in the studies
of tidal dissipation, and in the removal of tidal signals from their measurements prior to
any oceanographic and geophysical modeling (Shum et al., 1997).
In addition to practical applications such as navigation safety near the coasts and co-
tidal chart generation for mariners (Fang et al., 2004), ocean tide models also have a
variety of applications in scientific research. Examples include, but are not limited to, (a)

2
These are the comprehensive (i.e., both pelagic and coastal) tide gauge datasets for the entire world. The
locations of the tide gauges are shown in Figure 5.1.
2
study of regional tidal dynamics and dissipation (Han, 2000; Zu et al., 2008); (b)
operational ocean circulation modeling (Han et al., 2010), and climate forecast (Escudier
and Fellous, 2009); (c) inference of ocean tidal loading serving as a correction to GPS
and gravimetric measurement for study of regional dynamics of the solid Earth (Inazu et
al., 2009; Yuan et al., 2009) and its internal structure (Ito and Simons, 2011).
The accurate prediction of ocean tides serves as a critical correction for spaceborne
measurements. It is either directly provided as a geometric correction for the Earth-based
satellite altimetry data, or indirectly modeled as orbital perturbation sources for the
spaceborne gravity sensors, such as Gravity Recovery And Climate Experiment
(GRACE) or Gravity field and steady-state Ocean Circulation Explorer (GOCE) (Bosch
et al., 2009). These model corrections also enable an improved quantification and
modeling of global general ocean circulation (Fu and Cazenave, 2001), including
climate-sensitive signals of mass variations or transport. Most of the above-mentioned
applications require a gridded global ocean tidal forward prediction model.
While the ocean tide models derived from satellite altimetry have an accuracy of 2–3
cm RMS (root-mean-square) error in the deep ocean, their uncertainties inflate
significantly near coastal regions or over shallow seas and underneath ice-covered polar
ocean (Shum et al., 1997, 2001; Ray et al., 2011). The model accuracies are also region-
dependent, and vary among different models. A similar conclusion was found by Fok et
al. (2010).
There are two reasons for the above claims. Firstly, the temporal sampling rate of
altimetry data is insufficient which is limited by the satellite orbital design. The sampling
period at the same location is longer than the half daily and daily period of major ocean
tidal signals. This refers to tidal aliasing effect. Furthermore, when the data are sampled
very close to coasts, the data are further flagged by the contamination of the land surface,
and/or by the fact that the tide is on the ebb sometimes. As a consequence, this magnifies
the tidal aliasing effect, which is unfavorable for the accurate determination of ocean
tides using satellite altimetry. Secondly, the tidal prediction is complicated by the non-
linearity in dynamic tidal motion and local hydrodynamic processes. This is due to the
shape and bathymetry of the coast, in particular, when a complicated estuary system is
present. The overtide and compound tide are but one of the consequences that lower the
accuracy of ocean tides determination using satellite altimetry. The measurements from
coastal tide gauges are also prone to localized tidal effects which may not be
representative a few kilometers away from the coasts. In view of the above
considerations, ongoing investigation is necessary for a better data-constrained ocean tide
model with improved spatial and temporal resolution from multi-satellite altimetry,
particularly in the coastal regions.

1.2 Contemporary Ocean Tide Models Using Satellite Altimetry

The global ocean tide models can be classified into three groups: (1) hydrodynamic
model; (2) assimilation model; and (3) empirical model.
Hydrodynamic models are derived by solving the Laplace Tidal Equations (LTE)
using a finite difference or finite element method based on the bathymetry data around
3
the world, with the tidal heights from tidal constants of nearly all tide gauge records as
boundary conditions (or data constraints). Schwiderski’s (1980) model and the Grenoble
model derived by Le Provost et al. (1994) belong to this model type.
Assimilation models are derived either by assimilating altimetry data only, or in
combination with tide-gauge data, into the hydrodynamic model. The NAO.99b model
(Matsumoto et al., 2000) assimilates altimetry data only, whereas the TPXO6.2/7.1/7.2
(Egbert and Erofeeva, 2002) and FES94.1/99/2004 (Lefevre et al., 2000) models
assimilate both the altimetry data and the tide gauge data.
Empirical models are derived from altimetry data only using (a) response analysis
(Munk and Cartwright, 1966) extended with orthotide formalism (Groves and Reynolds,
1975; Cartwright and Ray, 1990b), (b) harmonic analysis, or (c) Proudman functions
(Sanchez and Pavlis, 1995). The empirical models could be subdivided into (i) semi-
empirical and (ii) purely empirical. The semi-empirical models use a-priori tidal
constants from a background ocean tide model, such as the FES model, for residual tidal
analysis (hereinafter called incremental tidal analysis) 3 . Examples for these models
include GOT00.2/4.7 (Ray, 1999), EOT08a/10a/11a (Savcenko and Bosch, 2008), and
DTU10 (Cheng and Andersen, 2011). The models derived from Cartwight and Ray
(1990a, b), Ray et al. (1994), Desai and Wahr (1995), and Smith (1999) are purely
empirical.
The hydrodynamic modeling approach provides robust information of the physics of
ocean tides and their dynamic characteristics in space. This approach includes the insight
and the understanding of tidal regimes, and the dependency on specific parameters such
as bottom topography, advection and bottom friction, tidal loading, and self-attraction.
However, their solutions are merely in agreement with in-situ observations in a
qualitative manner (Le Provost, 2001), because the approach heavily depends on the
availability and the quality of the bathymetry and tide gauge data for different regions.
Some fine (empirical) tunings of effects – such as a topographic drag and a bottom
friction coefficient – in specific regions are required to obtain a better solution
(Matsumoto et al., 2000; Arbic et al., 2004), and the cost of computation is enormous.
The assimilation modeling approach also suffers from the same weaknesses, as this
modeling approach exhibits a certain similarity with the hydrodynamic one, except for
the inclusion of altimetry data.
On the contrary, the empirical modeling approach offers less physical insight and
understanding into the tidal dynamics, but a simple, effective, and practical method for
accurate tidal analysis and prediction at a fixed location, provided that adequate altimetry
data observations are available. In addition, this approach also does not require tide gauge
data so that the accuracy comparison against tide gauges is independent. The accuracy of
ocean tide models computed from this approach is, in general, higher than that of the
hydrodynamic modeling approach.

3
Residual tidal analysis refers to solve for incremental tidal constants which is then added back to the
background tide model for the full ocean tide model. This is achieved by bilinear interpolating to the actual
locations of the satellite altimeter ground tracks from the nearby grids of the background tide model. Note
that the altimetry data used in the background tide model is also repeatly used in the residual ocean tide
modeling.
4
 Empirical models, no matter if purely empirical or semi-empirical, are derived using
two methods, in which a spatial smoothing step is required after the solution has been
made (Eanes and Bettadpur, 1995; Ray, 1999). The first method is to conduct a tidal
analysis at each location along the ground track of TOPEX satellite altimeter, followed
by a spatial interpolation onto a regular grid (Andersen, 1994, 1995a; Desai and Wahr,
1995) with homogeneous weighting applied to the corresponding altimeter data, which is
referred to as two-step method. When GFO and Envisat data are included, an iterative
solution step has to be made. This is achieved by using a predetermined TOPEX-alone
regular gridded tidal analysis solution as a background ocean tide model for the
incremental tidal analysis of GFO and Envisat altimeter data, because their tidal aliasing
effects are much worse than that of TOPEX. The tidal aliasing effects of GFO and
Envisat worse than that of TOPEX are due to their respective temporal sampling period at
17-day and 35-day interval, as will be illustrated in section 3.2 of this dissertation.
The second method is to acquire multi-satellite altimeter data at a search distance
from a predefined regular grid center location. The data are not reduced to the grid center
location a priori but are reduced in the least-squares solution process through weighting
the data in space via the Gaussian distance decay function (Smith, 1999; Savcenko and
Bosch, 2008; Bosch et al., 2009), which is referred to as one-step method. The underlying
reason of this method is to mitigate the tidal aliasing effect by the usage of more data
from multi-satellite altimeter with distinct spatial and temporal coverage. Most
contemporary models (i.e. GOT00.2/4.7, EOT08a/10a/11a, and DTU10) were made via
incremental tidal analysis based on FES tide model as background model.
Note that in the modeling process, the radial displacement of ocean tidal loading are
computed from direct computation of ocean tide estimates (Desai, 1996; Ray, 1999) in an
iterative fashion, from the 7% rule of ocean tides (Andersen, 1999) for the empirical
models, or they are computed from other ocean tidal loading models to correct the
altimetry measurements before empirical modeling.

1.3 Problem Statement and Research Methodology

In summary, the accuracy of 2–3 cm RMS (root-mean-square) error in the deep ocean
can be achieved for most ocean tide models, but notably larger uncertainties are found
among the ocean tide models near the coasts in the presence of a high hydrodynamic
and/or in the seasonally frozen (i.e. ice-covered) polar ocean. The error of the global tide
models against coastal tide gauges can exceed tens of centimeters or even a few meters at
some fixed locations. This is attributable to flagged altimetry data very close to the coast,
and the presence of non-linear dynamic tidal motion caused by local hydrodynamic
processes. The model accuracies are also region-dependent, and vary among different
models.
Even though the assimilation model, FES2004, has assimilated altimetry and 670 tide
gauges (all over the world) into the hydrodynamic model, the model accuracy near
coastal regions is still region-dependent even when compared with the same tide gauge
data used.

5
 The semi-empirical approach generates a full ocean tide models which are dependent
on tide gauge data, since FES model is mostly used as background ocean tide model for
incremental tidal analysis. As a consequence, this causes the final accuracy assessment
invalid to a certain extent when compared with the tide-gauge-derived harmonic
constants. In contrast, the purely empirical approach is a viable alternative because it
utilizes altimetry data only. This also allows an independent comparison of the derived
model with the tide gauge data.
This study aims to generate a purely empirical global ocean tide model with
0.25o×0.25o spatial resolution. The resulting solution is referred to as OSU12 ocean tide
model. A one-step orthotide response analysis is employed for the OSU12 model
generation via a novel spatio-temporal combination approach for data weighting and a
robust estimation technique for potential outlier downweighting, using improved multi-
satellite altimetry data from TOPEX, Jason-1/-2, Envisat, and GFO.
It is impossible to tell which process – spatial, temporal, or a combination of both – is
dominant in particular sea regions for modeling the remaining errors (i.e. unmodeled
ocean circulation signal, error in geophysical corrections and random noise, etc) that
behave as if stochastic noise. In contrast to previous one-step solution method for
empirical ocean tide modeling, a spatio-temporal combination approach is developed in
this dissertation that simultaneously considers both spatial and temporal (co-)variances
for weighting altimetry data in a balanced sense. Particular emphasis is paid on spatial
(co-)variances based on bathymetry depth for a substantial improvement along the
world’s coastal regions and over continental shelves when compared to existing models.
Besides the usage of a longer data time span (similar to current EOT08a/10a/11a and
DTU10 models, except Jason-2 data was excluded), and an improved quality of multi-
mission satellite radar altimetry (i.e., TOPEX/ Jason-1 and their tandem mission, GFO
and Envisat mission) data over the past few years, the potential improvement lies in an
appropriate spatial and temporal (co-)variance specification for the novel spatio-temporal
combination approach: (1) Spatial (co-)variance is designed via the modification of
dynamic interpolation method (Andersen, 1999) that considers the rapid distance decay
property near coasts in Shallow Ocean. This is because tidal dynamic characteristics
changes significantly when tidal wave from the Deep Ocean enters into the continental
shelf and shallow water regions where highly varying bathymetry depths are present; (2)
Temporal (co-)variance is assigned by considering the noise level of each altimeter for
each 1-Hz along-track measurement together with the total error budget from different
error sources. In contrast, previous tide models were generated by weighting each
observation homogeneously (i.e. equally weighted) or in space based on Gaussian decay
function with respect to grid location of the tidal solution.
To examine the effectiveness of the spatio-temporal combination approach, a
sensitivity analysis is conducted for an accuracy comparison against the equally weighted
solution, the weighted solution based on spatial (co-)variances, and the weighted solution
based on temporal (co-)variances under the same settings using independent tidal
constant from coastal tide gauges. Evaluations are made to the resulting ocean tide
model, OSU12, through (i) an external comparison with other contemporary models
using independent tidal constants from pelagic and coastal tide gauges around the world,
6
and (ii) an internal comparison via a variance reduction test using all altimetry data,
including independent Jason-2 altimetry data.
Chapter 2 introduces the basic theory of ocean tide generation and analysis
techniques. Chapter 3 introduces the principle of satellite altimetry, its data
characteristics and associated corrections, and the tidal aliasing effect due to under
sampling in terms of observed time interval. Chapter 4 provides a review of ocean tides
modeling approaches and the developed weighting technique via spatio-temporal
combination approach for the ocean tides modeling. Chapter 5 presents the detailed
analyses and the associated results, and the evidences of seasonal tidal variation in
Subarctic Ocean. Chapter 6 concludes the work in this dissertation.

7
 Chapter 2: Theory of Ocean Tides and Tidal Analysis
In this chapter, we review the theory of the gravitational forces of the Sun and the
Moon acting on the Earth that cause periodic deformations of the non-rigid Earth,
including the ocean, which are known as tides. The mathematical development of these
forces can be developed into a series of harmonic waves based on the tide-generating
potential will be presented in section 2.2. We describe, in brief, the procedures to derive
the amplitudes of the sea surface elevations associated with the tidal harmonics from
Doodson (1921) and Cartwright and Taylor (1971), and Cartwright and Edden (1973).
These classical treatises assumed that the Earth is comprised of only ocean, and that there
is an equilibrium response of the oceans to the tide-generating forces. The assumption is
not entirely correct for true ocean tides representation, for instance, in regions with high
dynamic oceanic variability where tidal mixing, internal tides, and barotropic response
for non-tidal features are dominant. However, the theory based on equilibrium tides does
serve as an important reference system for tidal analysis (Pugh, 2004).

2.1 Tide-generating Potential and its Species

According to the Newtonian law of gravitational attraction, the gravitational force
acting on a mass particle of the Earth at an arbitrary point P, near and outside the surface
of the Earth, with mass M P , by the mass of the perturbing terrestrial body, M , is given
by:
GM P M
F  (2.1)
R P2
where G is the universal gravitational constant and RP is the distance between the Earth
surface of the arbitrary point P and the center of the other body, either the Moon or the
Sun. Thus, the gravitational potential at point P due to the Moon or the Sun is
GM
UP  (2.2)
RP
This is the conventional definition in geodesy. Note that an alternative convention
involving a negative sign in equation (2.1) is adopted in physics. For the positive
definition of the gravitational potential in geodetic convention, an increase in potential
can be thought as an increase of the sea level of the ocean surface.

8
 Figure 2.1 The geometric concept for tidal harmonic potential development

Applying the cosine law for the EPB as shown in Figure 2.2, the distance between
the Earth surface at the arbitrary point P and the center of the perturbing body can be
represented as
1
 R R2 2
R P  r 1  2 e cos   2e  (2.3)
 r r 
Hence, the potential U P can be expanded into a series of Legendre polynomials with
increasing powers of Re / r , namely
1

GM  Re Re2  2
U P  U P ( Re ,  ,  )   
1 2 cos   2
r  r r 
(2.4)
GM  Re 2
R R3

 1  P1 (cos  ) e
2
P2 (cos  )  e
3
P3 (cos  )  
r  r r r 
where cos  is a function of  and  as provided in equation (2.7). We use Pl (cos  ) to
represent the Legendre polynomial of degree l . The tide-generating potential on the
Earth’s surface at point P caused by the Moon can be expressed as
l
GM   Re 
U ( Re ,  ,  )     Pl (cos )
r l 2  r 
(2.5)

where R e is the mean radius of the Earth, r and  are the geocentric distance and zenith
distance of the planetary body, respectively. Notice that the terms for l  0 and l  1 are
absent in equation (2.5) since the focus here is the deformation due to tides. Because the
ratio Re / r , denoting approximately the sine of the parallax of the attracting body (i.e., the
Moon or the Sun), is a very small value (Taff, 1985), the fourth term and all higher order
terms are neglected. With the substitution of the relevant expression of P2 (cos ) , the
main term ( l  2 ) in equation (2.5) is thus given as
9
 3GMRe2 1
U 2 ( Re ,  ,  )  3
(cos 2  ) (2.6)
4r 3

Figure 2.2 Geometric relationship between geocentric longitude and latitude ( , ) for P, right ascension
and declination position (  ,  ) of the planetary body B, and the zenith distance  and hour angle H.

The planetary body is usually at a different orbital plane inclined to the equatorial
plane of the Earth, in which its position (  ,  ) is kept changing such that the direction of
the generated tidal bulge changes instantaneously with a different hour angle H at the
observer position ( , ). Using spherical trigonometric formulas, one relates the above
quantities to the zenith distance  (Figure 2.3), which is given by
cos  sin sin   cos cos cos H (2.7)
where
H  G     (2.8)
in which  G is the Greenwich sidereal time. Thus, Pl (cos ) can be expanded (e.g.
Hobson, 1965) into
l
(l  m)!
Pl (cos  )   ( 2   0 m ) Plm (sin  ) Plm (sin  ) cos( mH ) (2.9)
m 0 (l  m)!
By substituting equation (2.9) into equation (2.5),
GM  Re l l (l  m )!
U   ( )  (2   0m )
r l 2 r m 0 (l  m )!
Plm (sin  ) Plm (sin  ) cos( mH ) , (2.10)

10
where Plm (x) denotes the associated Legendre function of degree l and order m, and  0m
is the Kronecker delta symbol defined as
1, when m  0
 0m   (2.11)
0, when m  0
Focusing on all the terms for degree l  2 with the substitution of the relevant
expressions for P20 (sin ) , P21 (sin  ) , and P22 (sin  ) , one gets the corresponding main
part of the tide-generating potential of the planetary body as
2
U 2  U 2 m
m 0
(2.12)
3GMRe2  1 
 (1  3 sin  )(  sin  )  sin 2 sin 2 cos H  cos  cos  cos 2 H 
2 2 2 2

4r 3  3 
The above derivation explicitly illustrates the physical reasoning for the differences in the
tidal species. One can see that the first term with m  0 is the long-period species, having
a period of half of year for the sun and half a month in case of the Moon since the
variations with the square of the sine of the tide-generating body’s declination; the
second term with m  1 is called the diurnal species, containing H and thus showing once
per day variation, while the third term with m  2 is called the semi-diurnal species,
containing 2H and thus showing twice per day variation.

2.2 Harmonic Expansion of the Tide-generating Potential

The tide-generating potential, as manifested from equation (2.12), is a function of the
time-dependent position of the attracting bodies (i.e., either the Sun or the Moon) that are
treated as point masses or bodies with uniform densities. The motion of the Sun and the
Moon can be described by six astronomical angles approximately linear with time, such
that the tide-generating potential can be developed into a series of corresponding
harmonics (Smith, 1999).
After the first development of harmonic expansion by Darwin (1883), Doodson
(1921) derived the first comprehensive expansion of the tide-generating potential, using
Brown’s lunar theory and with the position of the sun and moon with respect to ecliptic,
i.e. the plane of the Sun’s apparent orbit around the Earth, in terms of ecliptic longitude
and latitude. The expansion includes over 400 terms when compared to that with 39
terms as derived by Darwin (1883). It is currently expanded over 20000 terms
(Kudryavtsev, 2004). Each species U 2m in the 2nd degree term U 2 of the tide-generating
potential can be written as
U 2 m  G m ( )   k cos(  k (t )   k  m ) (2.13)
k

where Gm ( ) are Doodson’s geodetic coefficients with the first three terms (in case of
the Moon) given as

11
 3GMRe2 1 3 2
G0 ( )  (  sin  )
4R 3 2 2
3GMRe2
G1 ( )  sin 2 (2.14)
4R 3
3GMRe2
G2 ( )  3
cos 2 
4R
R is the mean distance from the moon (or the sun) to the Earth,  k are the coefficients of
the harmonic expansion, and  k are additive phase corrections, which are multiples of
 2 introduced to obtain a series of all positive coefficient k and cosine functions only
for attaining uniformity for a particular tidal constituent k (Casotto, 1989).
The Doodson arguments at Greenwich, k , are expressed as
 k (t )  A  Bs  Ch  Dp  EN ' Fps (2.15)
where N    N ,   t  s  h and t , s, h, p, N , p s are the fundamental angles that
represent Greenwich mean solar time, mean longitudes of the moon, the sun, the lunar
perigee, the lunar node, and the solar perigee, respectively. The variation of s, h, p, N and
p s can be expressed through polynomials of time T (in units of the Julian century) from
noon of December 31 1899 (UT), based on Brown’s lunar theory and Newcomb’s theory
of the sun (Doodson, 1921; Casotto, 1989; Smith, 1999), given as:
 tsh
s  270 .43659  481267 .89057T  0 .00198T 2  0 .000002T 3
h  279 .69660  36000 .76892T  0 .00030T 2
(2.16)
p  334 .32956  4069 .03403T  0 .01032T 2  0 .000010T 3
N  259 .18328  1934 .14201T  0 .00208T 2  0 .000002T 3
p s  281.22083  1.71902T  0 .00045T 2  0 .000003T 3
The relationship between the angles is illustrated in Figure 2.4.

12
 Figure 2.3 Relationship between Greenwich mean solar time, t , Greenwich mean lunar time,  , the
Moon’s mean longitude, s , and the Sun’s mean longitude, h (Smith, 1999).

The argument number ABC.DEF in equation (2.15) can be expressed in the combination
of
k1k 2 k3 .k 4 k5 k6  A( B  5)(C  5).(D  5)(E  5)(F  5) (2.17)
which is known as the Doodson number that denotes each tide constituent, such that the
Doodson arguments at Greenwich can be written in
 (t  t )
 k (t )   k (t 0 )   (2.18)
k 0

where the frequency of each tide constituent k,   , and the astronomical argument at
k

epoch t0 ,  k (t0 ) are given by:

  k   (k  5) s  (k  5)h  (k  5) p  (k  5) N   (k  5) p (2.19)
k 1 2 3 4 5 6 s

 k (t 0 )  k1 (t 0 )  (k 2  5) s(t 0 )  (k 3  5)h(t 0 )  (k 4  5) p(t 0 )

(2.20)
 (k 5  5) N ' (t 0 )  (k 6  5) p s (t 0 )
The Doodson number lays out a foundation to expressing astronomical arguments for
practical harmonic tidal analysis, described in section 2.3.2.

13
 Because of the availability of accurate lunar and solar ephemerides via numerical
integration of orbits of planetary bodies, the astronomical constants, and the advancement
in computers, the tide-generating potential has been recalculated with their expansion
(Cartwright and Tayler, 1971, and Cartwright and Edden, 1973) in the following form,
which follows the tidal response formulation (Munk and Cartwright, 1966) as follows:
  l * 
U  Re g  clm (t )Wlm ( ,  ) (2.21)
 l 2 m0 
where g is the Earth’s mean gravity acceleration on the surface and Wlm ( ,  ) are the
complex spherical harmonics defined as
2l  1 (l  m)!
Wlm ( ,  )  (1) m Plm (sin  )e im (2.22)
4 (l  m)!
*
The time varying complex conjugate coefficients c lm correspond to the Greenwich
equilibrium tide of degree l and order m . For tides, with l  2 , the coefficients are
c 2*m  (1) m  0 m  Bk e i (  k ( t )   k ) (2.23)
k

where B k represents the equilibrium tidal amplitude as tabulated in Cartwright and

m 
Taylor (1971), and Cartwright and Edden (1973). The term (1) 0 m is introduced to
obtain all positive amplitudes Bk and cosine arguments in the harmonic expansion. A
more details on the review and harmonic expansion development can be found in Wenzel
(1997).

2.3 Tidal Analysis

The above two sections illustrate the tide-generating potential as well as its harmonic
expansion. However, the equilibrium response of ocean tides may not be entirely correct
for true representation of ocean tides near coasts, although long-period tides having a
period of a month or longer should be expected to closely follow the equilibrium theory
(Lambeck, 1988).
In 1775, Laplace formulated a dynamic tidal theory, called Laplace tidal equations
(LTE), to describe the horizontal currents and associated tidal elevations of the ocean
caused by the tide-generating forces. These equations provide a physical theory to the
expression of the sum of both the solid-Earth and the ocean tides, which will be described
in detail in the following sections. Two practical methods for ocean tidal analysis, called
harmonic and response tidal analyses, respectively, will be described.

2.3.1 Laplace Tidal Equations

The Laplace tidal equations (LTE) describe the horizontal currents and associated
tidal elevations of the ocean as a result of the tidal forcing. The three equations, based on
the conservation of momentum (horizontal) and the conservation of mass (vertical), are
as follows (Laplace, 1775):
14
 u 1 
 2 e v sin    g c  
t Re cos  
v 1 
 2 e u sin    g  c    (2.24)
t R e 
 

1
 H b u    H b v cos    0
t Re cos     
The first two equations are the horizontal momentum equations that give the horizontal
acceleration of the water in the ocean as a consequence of the total tide-generating
potential  (Hendershott, 1981). The third equation is the mass conservation equation,
such that the net flux into or out of the area of local bathymetry depth H b is balanced by
a corresponding change in water level (e.g. Pugh, 1987). In the equations (2.24), e is the
Earth’s mean rotation rate in which the Coriolis parameter 2 e sin  is formed; u( , ,t )
and v( , ,t ) are the east and north component of the depth-mean velocity (i.e., the
horizontal tidal currents), and  c is the tidal height which is composed of ocean tides,
ocean tidal loading, and the tidal deformation of the solid earth. Note that the geocentric
tidal height is the sum of ocean tide and ocean tidal loading which is observed by the
satellite altimeter. The above equations are numerically solved for (u, v,  c ) driven by the
total tide-generating potential  . The inputs are (i) the bathymetry depth with respect to
the mean sea surface and (ii) the tidal height computed from harmonic constants derived
from tide gauge measurements that serve as open boundary conditions (or data
constraints); the output includes the ocean tides and the ocean tidal loading. The solid-
Earth body tide is relatively well known and therefore given.
To illustrate the concepts of all tidal components generated from the total tide-
generating potential, one should learn that the composition of the total tide-generating
potential  , truncated at degree 2, can be expressed as (Ray, 1998):
  U 2  U 2  U   U  (2.25)
where U 2 is the tide-generating potential, U 2 is the additional potential caused by the
redistribution of mass inside the Earth in response to the tide-generating potential, U  is
the potential induced by ocean tide, and U is an additional potential generated from the
ocean tidal loading, called the loading potential. Each component of  is described in
detail below for a clear explanation of the term g c   in equation (2.24).
With tide-generating potential U 2 , the tidal deformation of the solid earth given by
hU
e  2 2 (2.26)
g
The additional potential U 2 in response to the tide-generating potential can be
expressed as
U 2  k 2U 2 (2.27)

15
where h2 and k2 are the second-degree Love numbers, assumed to be independent of
frequency (Munk and MacDonald, 1975).
U  , being the potential induced by ocean tide, can be expressed in terms of spherical
harmonics as
U   g   l l (2.28)
l

3 w
where  l is the l-th degree spherical harmonic of the ocean tide with  l  , and
2l  1  e
 w and  e are the mean density of the ocean water and the Earth, respectively. This
potential causes the deformation of the solid earth in form of the ocean tidal loading:
 ol   hl' l  l , (2.29)
l

and U is an additional potential caused by the ocean load tides that can be expressed as
U   g  k l' l  l (2.30)
l

where hl' and k l' are the load Love numbers of degree l. Hence, the term g c   in
equation (2.24) is rewritten as:

g c    g  g e  U 2  U 2   g ol  U   U  
2 2 2 
 g  1  k  h U  g 1  k '  h '   .
l
l l  l l
(2.31)

Notice that, as the term g ol  U   U  in equation (2.31) depends on the solution of

the ocean tide, equation (2.24) has to be solved iteratively. Due to slow convergence, an
alternative procedure was proposed to approximate the term  1  k l'  hl'  l  l in
l
equation (2.31) by a simple scaling relationship to the tidal height. Accad and Pekeris
(1978) and Schwiderski (1983) suggested the relationship to be 0.085  or 0.07  ,
respectively. Note also that the non-linearity term such as the bottom friction and
advection terms have to be taken into account in equation (2.24) in the case of shallow
water.

2.3.2 Harmonic Analysis

Tides are periodic in nature, and it can be expressed as the sum of sinusoids. The
standard approach to determine the tides is referred to as ‘harmonic analysis’. At time t,
the harmonic expression of the ocean tidal height,  , at a location (φ, λ), can be
expressed as
 ( ,  , t )   H k ( ,  ) cos[  k (t )   k  G k ( ,  )] (2.32)
k

16
where k (t ) and  k share the same definition as in equations (2.13) and (2.23), and
H k ( ,  ) and Gk ( ,  ) are the unknown amplitudes and Greenwich phase lags of the
tidal constituent k which will be determined through a least-squares process.
Equation (2.32) is not in linear form. For a convenient least-squares estimation and to
avoid the singularities at the point of zero amplitude of any harmonic constituent of the
tide, which refers to as an ‘amphidromic point’, it can be decomposed into the cosine and
sine functions in linear form through the trigonometric addition formula as
 ( ,  , t )   [C k ( ,  ) cos(  k (t )   k )  S k ( ,  ) sin( k (t )   k )] (2.33)
k
where
C k  H k cos(Gk )
(2.34)
S k  H k sin(Gk )
which are called in-phase and quadrature amplitude terms, respectively, that constitutes
the tidal (harmonic) constants 4 . Their relation to the amplitude and phase lag is given by
H k  Ck2  S k2
Sk (2.35)
Gk  arctan( )
Ck
Because the lunar nodal tidal effect occurs with a cycle of 18.61 years, the nodal
correction has to be applied to take into account of the slow variation for the longitude of
the lunar node. Notice that it is strictly to the lunar tidal constituents. Therefore, the lunar
nodal factor, f k , and the nodal angle,  k , are introduced to the harmonic expression of
the tidal height in equation (2.33) (e.g., Munk and Cartwright, 1966; Schureman, 1971)
and expressed as
 ( ,  , t )   f k C k ( ,  ) cos(  k (t )   k   k )  f k S k ( ,  ) sin( k (t )   k   k ) (2.36)
k
Both of them depend on the longitude position of the lunar node with slow variation in
time throughout the cycle.

2.3.3 Response Analysis

Notwithstanding the convenience of harmonic analysis for the harmonic constants
(i.e., C k , S k ), it has the inability of resolving frequencies close to each other due to
limitation of the data time span. Munk and Cartwright (1966) developed the response
analysis method that relates the equilibrium tide, c2*m (t ) , (an input), to the ocean tidal
height at a fixed location,  ( ,  , t ) , (an output), by the response weight function ,
w2 m ( ,  , t ) , (a system), with the assumption of the credo of smoothness in which sharp
resonance peaks are presumably not occurred for the ocean responses to gravitational
4
Note that the terminologies, (i) amplitude and phase and (ii) tidal (harmonic) constants refer to the same
quantities which are always interchangeable in this context.
17
forcing. This idea, which is originated from electrical engineering, has a conceptual
advantage that distinguishes the equilibrium tide from the ocean dynamic system. It also
permits an analysis in spectral domain. Though it is fundamentally designed for linear
systems, it can be extended to non-linear equations that describe shallow water tide
propagation (Pugh, 1987). Therefore, the ocean tidal height is expressed as a convolution
of the equilibrium tide, c2*m (t ) , and the response weight function, w2 m ( ,  , t ) .
Considering only the main tides ( l  2 ), the ocean tide height  (the output) at time t
and at location (φ, λ) is given by
 2 * 
 ( ,  , t )  Re c2 m (t ) * w2m ( ,  , t ) (2.37)
m0 
where
c 2* m (t )  a 2 m (t )  ib2 m (t ) (2.38)
is the equilibrium tide of second order, expressed in an equivalent form to equation
(2.23), where a2m (t ) and b2 m (t ) are the real and imaginary parts of the equilibrium tide
representation that contain the nodal parameters implicitly and, hence, no nodal
correction factor has to be included in the analysis. The response weight function (i.e., the
system) can also be represented as
S
w2 m ( ,  , t )  w
s  S
2 ms ( ,  ) (t  sT ) (2.39)

where  (t ) is the unit impulse function. T is the time lag interval, usually taken as two
days. S is the number of response weights which is generally be chosen as 1 or 2. Thus,
the response method represents the ocean tides as a weighted sum of past and future
values of the equilibrium tide. Note that the use of future values of equilibrium tide
(negative values of s) is reserved for a better illustration of the mathematical concept.
The response weight w2 ms can be further written as
w2ms ( ,  )  u2ms (,  )  iv2ms ( ,  ) (2.40)
such that equation (2.37) can be rewritten as
2 S
 ( ,  , t )    [u 2 ms ( ,  ) a 2 m (t  sT )  v 2 ms ( ,  )b2 m (t  sT )] (2.41)
m 0 s  S

where u2ms and v2ms are the response weights that will be determined from the ocean tidal
height observations through a least-squares process. Note that the response weights
define the admittance function of each tidal constituent with angular frequency  
k

through a Fourier transform of itself, which is given as

 )  X (  )  iY (  )   w ( ,  , t )e i k t dt
Z (
2m k 2m k 2m k 

2m (2.42)
such that the response weights are represented in the spectral domain. Inserting equations
(2.39) and (2.40) into equation (2.42), one obtains
S
 )  (u

Z 2 m ( k 2 ms  iv 2 ms )e i k sT (2.43)
s  S

18
with
S
 )
X 2 m ( k [u
s  S
2 ms
 sT )  v sin(
cos( k 2 ms
 sT )]
k

S
(2.44)
 )
Y2 m ( k  [v
s  S
2 ms
 sT )  u sin(
cos( k 2 ms
 sT )]
k

It should be noted that the time lag interval T and the number s govern the degree of
smoothness which should only be applied to luni-solar gravitational tidal forcing on the
ocean (i.e., dominant tides) according to Munk and Cartwright (1966). The obtained
admittance function allows a number of small tides to be inferred by interpolation in the
spectral domain which is assumed to be well defined by the dominant tides, thereby,
making it easier to separate tidal constituents with very close frequencies than that of the
harmonic analysis. Though it has fewer parameters to be solved when compared to
harmonic analysis, the stability of the tidal solution is, in general, better.
Equation (2.44) yields the admittance parameters that provide an equivalent
expression to describe the tides in the spectral domain as the harmonic coefficients would
in equation (2.34). The relation between admittance and harmonic coefficients is given by
(Cartwright and Ray, 1990b; Smith, 1999):
 )
C k  (1) m 0 m Bk X 2 m ( k
(2.45)
S  (1) m  0 m
B Y (  )
k k 2m k

where Bk is the equilibrium amplitude in (2.23).

Notice that the functions a 2 m and b2m in equation (2.38) are not orthogonal, which
may result in an ill-conditioned normal matrix when it is based on (2.41).

2.3.4 Orthotide Formulation and its Relation to Harmonic Constants

To ensure orthogonal basis functions and to avoid an ill-conditioned normal matrix,
the orthogonalized convolution method for the tidal prediction was introduced by Groves
and Reynolds (1975). This leads to the so-called orthotide formulation that exhibits
similarity with equation (2.41), in which the tidal height can be written as
2 2S
 ( , , t )  [U mj ( ,  ) Pjm (t )  V jm ( ,  )Q mj (t )] , (2.46)
m 0 j 0

where usually S  1 is adopted since it is shown to be sufficient in practice for accounting

most tidal variation at several sites (Alcock and Cartwright, 1978; Cartwright and Ray,
1990b). U mj ( ,  ) and V jm ( ,  ) are the unknown coefficients of the orthotide functions
to be determined via a weighted least-squares adjustment. Pjm (t ) and Q mj (t ) are the
orthotide functions, defined by linear combinations of a2 m (t  sT ) and b2m (t  sT ) ,
that generate the coefficients for the design matrix with the expressions for the first few
orthotide functions ( j=1,2,3), following Cartwright and Ray (1990b), represented as:

19
 P0m (t )  p 00
m
a 2 m (t )
P1m (t )  p10m a 2 m (t )  p11m a 2m (t )
m  m 
P2m (t )  p 20
m
a 2 m (t )  p 21 a 2 m (t )  q 21 b2 m (t )
(2.47)
Q0m (t )  p 00
m
b2 m (t )
Q1m (t )  p10m b2 m (t )  p11m b2m (t )
m  m 
Q2m (t )  p 20
m
b2 m (t )  p 21 b2 m (t )  q 21 a 2 m (t )
where the first few coefficients of pijm and q ijm are listed in Table 2.1.
Table 2.1 Orthotide coefficients

Diurnal Semidiurnal
(m=1) (m=2)
p00 0.0298 0.0200
p10 0.1408 0.0905
p11 -0.0805 -0.0638
p20 0.6002 0.3476
p 21 -0.3025 -0.1645
q 21 0.1517 0.0923

Here, a 2m (t ) , a 2m (t ) , b2m (t ) , and b2m


(t ) are defined as:
a2m (t )  a2m (t  T )  a2m (t  T )
(2.48)
b2m (t )  b2 m (t  T )  b2 m (t  T )
Expression for the transformation from U mj and V jm to the admittance is given by
Cartwright and Ray (1990b):
X 2 m (  )  p m U m  [ p m  p m c ( )]U m  [ p m  p m c( )  q m s ( )]U m
k 00 0 10 11 k 1 20 21 k 21 k 2
(2.49)
 m m 
Y ( )  p V  [ p  p c( )]V  [ p  p c( )  q s( )]V
m m m m m  m  m
2m k 00 0 10 11 k 1 20 21 k 21 k 2
where
 )  2 cos(
c (  T )
k k
(2.50)
 
s( )  2 sin( T )
k k

2.4 Chapter Summary

We briefly illustrated the fundamental forces of the terrestrial bodies (such as the Sun
and the Moon), exerting on the Earth and thereby generating tides. The forces can be
expanded into a series of harmonics based on the tide-generating potential. Laplace tidal
equation is introduced as a fundamental basis for physical understanding of tides and its
dynamics. For practical tidal estimation and prediction, harmonic analysis is described as

20
the standard method. We also introduce the method of response analysis for the response
weight that allows tide to be inferred via interpolation in the spectral domain which can
be later converted back to tidal constants. Its orthogonal representation is also described.
Due to the aforementioned disadvantage of inability to resolve tidal frequencies close
to each other in the harmonic analysis method, the response analysis method that ensures
the orthogonality of the basis functions and avoids any ill-conditioned normal matrix is
employed in this dissertation.

21
 Chapter 3: Satellite Altimetry Technique
Before the advent of satellite altimetry, ocean tides have historically been measured
along the coastlines by tide gauges and at islands either by tide gauges or by bottom
pressure recorders. Hence, the empirical knowledge of global ocean tides in the deep
ocean remained largely unknown.
The advent of satellite altimetry enables the global ocean tides to be observed in the
deep ocean. It also addresses a vast number of fundamental and interdisciplinary
scientific questions, thanks to the satellite orbital design allowing synoptic means of
global ocean surface (or land surface) height observations with an approximate weekly
(or longer) temporal sampling. The accuracy of satellite altimeter measurement has been
improved since the launch of TOPEX/POSEIDON in the 1990s, primarily due to
precision orbit determination techniques, gravity field and other force modeling, along
with the improvement in instrumental, propagation media (e.g., troposphere and
ionosphere), and geophysical corrections.
In this chapter, we illustrate the difference of measurements through satellite
altimetry when compared to those by tide gauges, particularly in terms of spatial and
temporal data sampling. We review a measurement principle of satellite altimetry,
together with corrections that have to be applied to the altimeter measurements, in section
3.1. Sampling aliasing effect, in particular tidal aliasing, is introduced to gain a better
understanding on the temporal sampling aspects of altimeter data. This will be discussed
in section 3.2.

3.1 Sea Surface Height measurement in satellite altimetry

Before the era of satellite altimetry, tide gauges have been the only data source for
ocean tide modeling. They are mostly located near coastal regions for the accurate
measurement of ocean tides and sea level, in which the measurements are ultimately for
maintaining the ship cargo safety and other navigation application purposes.
Sea level observations from tide gauges are taken at a regular interval of 5, 10, 15,
and 30 minutes. This temporal sampling rate of sea level allows accurate determination of
ocean tides. However, due to their locations along coasts and the sparseness of the global
tide gauge network, tide gauge measurements provide limited knowledge about the
spatial characteristics of tides in the open ocean.
Compared to tide gauges, satellite altimetry measures sea surface height for
monitoring global sea level variation and ocean tides in a synoptic manner. However, it
has an approximate weekly (or longer) temporal sampling rate at certain fixed locations
on the Earth, due to the orbital design of the so-called Exact Repeat Orbits for satellite
altimetry.
The main difference of tide gauges and satellite altimetry lies in both the spatial and
the temporal sampling characteristics (Figure 3.1). It is apparent that the temporal
sampling of the tide gauge data is more than sufficient to separate the tides well from
other oceanographic signals of interest along with data noise (Parke et al., 1987), when

22
compared to that of satellite altimetry. However, sea level measurements are sampled
globally using satellite altimetry.
Satellite altimetry uses radar to measure the round trip travel time taken by a radar
pulse to the sea surface and return back to the satellite receiver. Because of the favorable
property of a relatively flat water surface with waves, the pulse-limited radar altimetry is
designed to be especially operated over the ocean. Measurement principles of satellite
altimetry are described in the following subsection.

23
 2

Figure 3.1 Tide gauge locations (red triangle) and theoretical altimeter tracks (blue – TP; green – GFO; red
- Envisat) [top panel] in Gulf of Mexico region and the time series of both tide gauge and TOPEX altimetry
measurements around the tide gauge, 1993–1999, in the Louisiana wetlands (yellow squares) [bottom
panel].

24
3.1.1 Measurement principle of satellite altimetry
The concept of satellite altimetry was first envisioned at NASA’s Williamstown
Conference in 1969 (Kaula, 1969) and the technique of satellite altimetry had,
subsequently, been demonstrated a few years later during the SKYLAB missions in
1970s (Seeber, 1993). Following the success of this experiment, new and improved
altimetry missions have been launched over the past two decades. Table 3.1 lists the basic
information about past, present, and planned future satellite radar and laser altimetry
missions, along with their orbital characteristics and repeat-orbit periods.

Figure 3.2 The schematics of satellite altimetry observations (Courtesy of AVISO)

25
 Table 3.1 Satellite radar and laser altimetry missions

Inclination Repeat Period

Mission Launch period
(degree) (days)
GEOS-3 Apr 1975 – Dec 1978 115 N/A
SEASAT Jul 1978 – Sept 1978 108 17.0505
SEASAT Sept 1978 – Oct 1978 108 3
GEOSAT GM Mar 1985 – Nov 1986 108 N/A
GEOSAT ERM Nov 1986 – Dec 1989 108 17.0505
ERS-1 July 1991 – Jun 1995 98.5 3, 35, 168 a
TOPEX/POSEIDON Aug 1992 – Dec 2005 66 9.9156
ERS-2 Apr 1995 – Jun 2003 98.5 35
GFO May 1998 – Sept 2008 108 17.0505
JASON-1* Nov 2001 – present 66 9.9156, 406 b
ENVISAT Mar 2002 – present 98.5 35
ICESat-1 Jan 2003 – Aug 2010 94 8, 91 c
JASON-2 June 2008 – present 66 9.9156
CryoSat-2 Apr 2010 – present 90, 92 c 2, 369 d
HY-2A Aug 2011 – present 99.35 14
SARAL /AltiKa e Sum 2012 98.55 35
Sentinel-3 2013 98.65 27
Jason-3 Apr 2014 66 9.9156
ICESat-2 Early 2016 92 91
a. The three-day repeat periods refer to the mission phases for the calibration and ice observations primarily
following the InSAR observational requirements, the nominal mission orbit, and the geodetic phase,
respectively.
b. JASON-1 has interleave/tandem/geodetic mission phases.
c. The validation orbit and the mission orbit have the 8-days and 91-day repeat period, respectively.
d. The validation orbit and the mission orbit have their inclination and repeat period of 92-deg and 369 days
with sub-cycle of 30 days.
e. Satellite with ARgos & ALtika (SARAL) is a French-Indian mission for the monitoring of the
environment: Altimetry (AltiKa) and contribution to ARGOS system (http://smsc.cnes.fr/SARAL/).

The measurement principle of an Earth orbiting satellite radar altimeter is

geometrically illustrated in Figure 3.2 which can be mathematically expressed as (Parke
et al., 1987; Wagner, 1989):
hssh  horb  halt  hMSS  hs (3.1)
where hssh is the sea surface height with respect to a predetermined reference ellipsoid,
horb is the satellite orbital altitude with respect to the International Terrestrial Reference
Frame (ITRF) and relative to the reference ellipsoid, halt is the altimeter range (or height)
from the satellite to the sea surface that is determined by multiplying the speed of light
with a half of the measured two-way travel time of the radar/laser pulse transmitted from
the altimeter antenna and reflected by the sea surface, hMSS is the mean sea surface
(MSS) comprising geoid undulation and dynamic sea surface topography, hs is the

26
(instantaneous) sea surface height (SSH) anomaly 5 . The orbit is determined by satellite
tracking data including Doppler Orbit determination and Radiopositioning by Satellite
(DORIS), Satellite Laser Ranging (SLR), Precise Range And Range-Rate Equipment
(PRARE), and Global Positioning System (GPS).

Figure 3.3 Sea surface height (SSH) anomaly residuals, after instrumental, geophysical, and media
corrections, were applied, including the ocean tide correction using NAO.99b (Matsumoto et al., 2000)
ocean tide model.

The main goal of satellite radar altimetry is to study the general ocean circulation and
to enrich the understanding of the role of ocean circulation in Earth’s climate at global
(i.e., gyre and basin) scales (Fu et al., 1994). Other applications of the altimetry
measurements include ocean tides, marine geodesy and geophysics, ocean wave height,
wind speed (Chelton et al., 2001), hydrology, ice-sheet elevation and sea-ice freeboard
elevation changes, and even solid-Earth geodynamics applications (Lee, 2008). When the
tidal fluctuations and sea level variations due to changes of solar heating, atmospheric
pressure, and wind are excluded, the precise instantaneous sea surface height, hssh ,
obtained from altimetry measurements represents the sum of MSS and the SSH anomaly
residual 6 (Calman, 1987).

5
In oceanography, Sea Surface Height (SSH) anomaly refers to the Mean Sea Surface subtracted from
SSH. In geodetic science, anomaly is usually thought of as gravity anomaly, which is the difference in
gravity on geoid at point P and ellipsoid at point Q.
6
Sea Surface Height (SSH) anomaly residuals, after all known geophysical corrections are applied, are the
general ocean circulation signals that are subject to further modeling. Residuals refer to the prediction of
the random errors after the least-squares adjustment solution process.
27
 Figure 3.3 illustrates the spatial variation of SSH anomaly residuals after
instrumental, geophysical, and media corrections were applied, along with the subtraction
of mean sea surface (MSS). Substantial ‘non-tidal’ dynamic oceanic variability is present
for some coastal regions with a well-known feature, called ‘west boundary current’, due
to the general ocean circulation.

3.1.2 Corrections to altimeter radar measurements

Direct altimetry measurements from an Earth-orbiting platform require several
corrections to achieve the desirable sea surface or land surface topography measurements.
Orbit error is one of the major errors resulting from inaccurate Earth’s gravity model used
for precise orbit determination (POD) (Tapley and Rosborough, 1985), even during or
before the TOPEX/Poseidon era. Other errors include (conservative and non-
conservative) forcing models (e.g. atmospheric drag, solar and Earth radiation pressure),
solid Earth and ocean tides, satellite-originated thermal forces, the measurement errors in
various tracking systems (e.g., DORIS, SLR, PRARE and GPS) such as measurement
sampling, atmosphere delay, and other related sources, and the errors in the terrestrial
reference frame and its realization. They all contribute to the total error budget.
The effect of any orbit error must be removed from the satellite altitude before it is
used to derive the sea surface height (SSH). Current POD accuracy for a typical modern
satellite altimeter, e.g., Jason-2, is better than ±2 cm (Bertiger et al., 2010), even for the
near-real-time (NRT) operational Geophysical Data Record (GDR) products (Desai and
Haines, 2010). The standard error of ±2 cm is adequately small for ocean tide modeling.
Hence, equation (3.1) is a good representation for hssh provided that the altimeter height
measurement, halt , has been corrected for systematic effects. One should also notice
other innovative orbit determination procedures, such as the reduced-dynamics approach
(Bertiger et al., 2010) and the purely geometric approach (Kwon et al., 2003), as opposed
to the dynamical approach (Tapley et al., 1994a,b). The reduced-dynamics and the
purely geometric approaches are based on the fact that there are abundant tracking data
available, i.e., high-low GPS to Low Earth Orbiters (LEO) or altimeter satellite tracking
data, and DORIS tracking data of modern altimetry satellites.
To obtain accurate (i.e., rather bias-free) SSH observations from satellite altimetry,
three major kinds of corrections have to be applied to the altimeter height measurements.
They are classified as: (i) instrument corrections, (ii) media propagation corrections, and
(iii) geophysical corrections. In addition, the instrument bias and the corrections are
calibrated using either absolute calibration sites (Haines et al., 1990; Christensen et al.,
1994; Shum et al., 2001), or relative calibrations are performed by comparing to other
concurrently flying satellite altimeters. These three types of correction are elaborated in
details below.

(i) Instrument corrections:

28
 Instrument corrections refer to sensor related corrections. They include Doppler-shift
correction, center-of-mass offset, nadir pointing error, instrument temperature and clock
(Ultra Stable Clock, or USO) related corrections, and some internal calibration biases.
Doppler-shift is due to the Doppler frequency shift, caused by the radial velocity of
the satellite. It affects the time delay measurement and, therefore, the range.
Center-of-mass offset accounts for the difference between the phase center of the
altimeter antenna, where the radar pulse is both transmitted and received, and the mass
center of the satellite on which the orbit computation is based. The satellite mass center
changes, depending on the spacecraft, as a function of fuel usage, and required operations
of a particular satellite including solar panel motions, and orientation changes (i.e.
satellite yaws) of the satellite.
Nadir pointing corrections refer to the offset caused by the deviation of the beam
direction from the vertical; the range measurement, hence, results in a slant-range to the
point offset from the nadir point.
Besides the aforementioned instrumental errors, the bias of the altimeter range
measurement toward the troughs of ocean waves, due to the sea state, should be taken
into account. It arises from three interrelated effects: tracker bias, skewness bias, and
electromagnetic (EM) bias (Rummel, 1993). Sometimes, the sea-state bias is categorized
as a geophysical correction.

(ii) Media corrections:

Media corrections refer to the media propagation correction due to radar pulse
passing through ionosphere and troposphere before reaching the sea surface (Figure 3.1).
They include the ionospheric correction and, both dry and wet tropospheric corrections.
The ionosphere delay is a frequency-dependent correction, in which the effect of the
first-order ionosphere correction is between 5 cm and 20 cm given the frequency domain
of 14GHz (Lorell et al., 1982). It is generally corrected by combining the onboard dual-
frequency altimeter measurements.
The dry tropospheric correction is a correction for the dry-air component in the
troposphere. As it cannot be measured directly by sensors onboard the satellites, the dry
troposphere delay is usually corrected by models such as the one by Saastamoinen
(1972).
The wet tropospheric correction is a correction for the integrated water vapor content
in the troposphere, which is comparatively more difficult to model than its dry-air
component counterpart. It can be corrected either using direct measurements from the
onboard microwave radiometer (Tapley et al., 1982), or from models such as the
European Centre for Medium-Range Weather Forecasts (ECMWF) model (Faugere et al.,
2011).

(iii) Geophysical corrections:

Geophysical corrections refer to systematic geophysical effects that can be modeled
and, consequently, be corrected. These consist of the barotropic response of the ocean to
the atmospheric pressure, i.e., an inverted barometer (IB) model, and corrections for

29
various tidal effects including the solid Earth body tides, ocean tides, ocean tidal loading,
and the pole tide.
The atmospheric barotropic correction is to remove the effect due to the ocean surface
deformation under the atmospheric loading. The inverted barometer (IB) model says that
1 millibar increase in atmospheric pressure will result in 1 cm decrease in the ocean
surface height (Ponte et al., 1991; Dorandeu and Le Traon, 1999).
The solid Earth body tide correction accounts for the periodic variations in both land
and sea surface, due to the tidal deformation of the underlying non-rigid Earth, including
the sea-floor, under the attraction of the astronomical bodies (Moon and Sun) (Munk and
MacDonald, 1975). It can be derived from the tide-generating potential, introduced in
Chapter 2 using Love numbers with the assumption of an elastic Earth with uniform
density of mass. Detailed information can be found in Chapter 6 of the IERS Conventions
(2010) and in the papers by Cartwright and Taylor (1971) and Cartwright and Edden
(1973). Incorporation of physically plausible assumptions about a lateral inhomogeneity
for the 3-D elastic Earth model has currently been proposed along with a viscoelastic
structure (Latychev et al., 2009); however, the result for the maximum change in the
amplitude under this assumption is below the 1 mm level which should be negligible in
view of the altimetry measurement quality of ±2 cm level.
Ocean tides represent a substantial time-variable component that is responsible for
sea surface deformation. The correction for ocean tides can be computed from some
available forward tide models, which will be introduced in section 4.1. Ocean tides also
cause an oceanic mass redistribution with the associated load change on the crust,
thereby, producing time-varying deformations of the earth. This is called the ocean tidal
loading effect (Ray, 1998). Similar to the solid Earth body tide, the displacement of the
earth crust caused by ocean tidal loading can be derived from the tide-generating
potential (Cartwright and Taylor, 1971; Cartwright and Edden, 1973). Schwiderski
(1980) proposed the 7% rule for ocean tidal loading, meaning that the ocean tidal loading
height is about 7% of the corresponding ocean tidal height, though a value of 8.5% was
suggested previously by Accad and Pekeris (1978). The ocean tidal loading correction,
first theoretically formulated by Farrell (1972), can be computed numerically from major
constituents of the ocean tides by using available program (Agnew, 1997).
The pole tide is a tidal response of both solid-Earth and the oceans to the centrifugal
potential caused by small perturbations to the Earth’s rotation axis (known as polar
motion). The perturbations occur at an annual period and 433-day period (also called the
Chandler wobbles) with amplitude of ~6 mm (Wahr, 1985; Desai, 2002). Its correction
can be computed if the location of the pole, as a function of the polar motion angles, is
known.

The observed sea surface height and its anomaly

With the introduction of the above corrections to the altimeter range measurements,
one can rewrite equation (3.1) in a more complete form as
obs
hssh : hssh  e  hMSS  hs  e
(3.2)
 hMSS  hb  ht  hinstr  hssb  hiono  hdry  hwet  hsol    hol  h pole  hIB  e

30
where hb is the potential satellite bias, ht is the height due to sea level trend, hinstr is the
instrument correction, hssb is the sea state bias correction, hiono is the ionosphere
correction, h dry is the dry troposphere correction, hwet is the wet troposphere
correction, h sol is the solid earth tide correction,  is the ocean tide correction (which is
the signal of interest and thus to be modeled in this dissertation), hol is the ocean tidal
loading correction, h pole is the pole tide correction, hIB is the inverted barometer
correction; e includes the error in the altimetry measurement, the error in geophysical
corrections, and the (unmodeled) general ocean circulation signal. Although the use of
existing numerical general ocean circulation model is a possibility for the removal of
general ocean circulation signal, it is subject to further scientific investigation before
serving as a geophysical correction because altimeter data has already been assimilated
into the general ocean circulation model that makes itself no longer independent of
altimeter data. Desai (1996) conducted an experiment to analyse the effect of including
general ocean circulation model (without altimetry data assimilation) as a correction to
altimeter measurement for empirical ocean tide modeling, but the effect is negligible. A
similar conclusion was found when an experiment is conducted using ECCO ocean
circulation model (KF080) in the Gulf of Mexico region. These findings are consistent
with other research contributions; indicating that the statistical distribution of the SSH
anomaly residuals behaves similar to stochastic noise (Desai, 1996; Niedzieski and
Kosek, 2009). Therefore, the error in geophysical corrections, the unmodeled general
ocean circulation signal and measurement noise are regarded as stochastic noise for the
ocean tide modeling.
In satellite altimetry context, the media and geophysical corrections can be served as
signal to be modeled or to be removed, depending on a particular research interest. In this
dissertation, since the ocean tide effect is the signal of interest, it should be regarded as
signal in the observations and not be treated as a correction.
Based on equation (3.1) and (3.2), the observed instantaneous sea surface height
obs
(SSH) anomaly, hssha , from satellite altimetry, after applying the instrument correction
and the environmental corrections except that for the ocean tides, can be expressed as
obs
h ssha  h ssh
obs
 h MSS  hinstr  h ssb  hiono  hdry  h wet  h sol  hol  h pole  h IB
(3.3)
 hb  ht    e
Note that equation (3.3) should be applied to the same along satellite track location
(  ,  ) only. In this dissertation, our solution approach used the original hssha obs
at (  ,  ) at
time epoch t to estimate the ocean tides at a predefined grid center location (  c , c ).
Hence, one should take into account both the (co-)variance in space and in time
simultaneously for data weighting using one-step empirical method, which will be later
described in section 4.2.

31
3.2 Tidal Aliasing
The estimation of ocean tides using satellite altimetry is limited by the temporal
sampling property due to the fundamental satellite orbit design. Because all altimetry
satellites have a repeat periods of a few days or more (i.e., 9.9156-day for TOPEX/Jason-
1/Jason-2; 17.0505-day for GFO; 35-day for Envisat), the diurnal (daily) and semi-
diurnal (half-daily) tides at a fixed location appears as long-period signals. This refers to
tidal aliasing effect.
Shannon sampling theorem stated that it is necessary to sample the signal at a rate at
least twice the frequency f k of the signal, i.e. f N  2 f k , where f N is the Nyquist
frequency, in order to fully reconstruct the original analog signal. In other words, a signal
with period Tk can be fully reconstructed if the sampling values are obtained at an
interval of less than Tk 2 , at least. Otherwise, the signal of period Tk will be aliased to a
longer period  P , which is referred to as aliased period. As a consequence of the aliasing
effect, one has to wait  P instead of Tk to fully reconstruct the signal of frequency f k .
Since the presence of data noise, unmodeled ocean circulation signal, and potential
data gaps due to flagged data for particular cycles in altimetry data, the observations have
to be collected more than just an aliased period  P .
The tidal aliased period for satellite altimetry with repeat period P can be computed
as below (Parke et al., 1987). The phase change, P , of a given tidal constituent k of
period Tk within the range   ,   , is given by
2  P
 P  (3.4)
Tk
such that the resulting principal alias period,  P , for a given tidal constituent k is
2  P
P  (3.5)
 P
An alternative formula for the aliased period can be found in Wang (2004).
Apart from the aliasing effect, another issue is the separation of different tidal
constituents with very close frequencies. The criterion, known as Rayleigh criterion, for
separating two neighboring tides is that they should at least differ in phase by a cycle
with a minimum required data time span, Tr (also known as Rayleigh period). It is given
by (Smith, 1999):
Tr 1  2  2 (3.6)
where 1 and 2 are the aliased angular frequencies. Thus, it can be simplified as
1 1 1
  (3.7)
Tr T1 T2
where T1  2 1 , T2  2  2 .

32
Table 3.2 Tidal aliased periods (diagonal) and Rayleigh periods (off-diagonal) in days for TOPEX/Jason-
1/-2, GFO, and Envisat, with repeat periods P of 9.9156, 17.0505, and 35 days, respectively. (Courtesy of
Smith, 1999)

M2 S2 N2 K2 K1 O1 P1 Q1 Mf Mm Ssa Sa
TOPEX/Jason-1/-2
M2 62 1084 245 220 97 173 206 594 87 50 94 75
S2 59 316 183 89 206 173 384 94 52 87 70
N2 50 116 69 594 112 173 134 62 68 57
K2 87 173 97 3355 349 62 40 165 114
K1 173 62 183 116 46 33 3355 329
O1 46 94 134 173 69 61 52
P1 89 316 61 40 173 118
Q1 69 76 46 112 86
Mf 36 116 45 40
Mm 28 33 30
Ssa 183 365
Sa 365
GFO
M2 317 361 62 121 393 175 341 97 88 52 431 2407
S2 169 75 183 4464 341 175 132 116 61 2232 314
N2 52 128 74 97 53 175 215 317 73 61
K2 88 175 393 90 475 317 91 169 115
K1 175 317 183 128 113 60 4464 338
O1 113 116 215 175 74 296 164
P1 4464 75 70 45 190 398
Q1 74 954 113 125 93
Mf 69 128 110 85
Mm 45 59 51
Ssa 183 365
Sa 365
Envisat
M2 95 95 3169 196 128 365 128 328 519 349 196 128
S2 ∞ 97 183 365 75 365 133 80 130 183 365
N2 97 209 133 328 133 365 446 393 209 133
K2 183 365 128 365 487 142 446 ∞ 365
K1 365 95 ∞ 209 102 201 365 ∞
O1 75 95 173 1236 179 128 95
P1 365 209 102 201 365 ∞
Q1 133 201 5253 487 209
Mf 80 209 142 102
Mm 130 446 201
Ssa 183 365
Sa 365

Table 3.2 illustrates the tidal aliased periods and Rayleigh periods of the most
prominent tidal constituents for TOPEX/Jason-1/-2, GFO, and Envisat, respectively. The

33
aliased periods are shown in the diagonal whereas the off‐diagonal elements are the

Rayleigh period indicating the minimum data time span required for separating the tidal
constituents from each other. The tidal signal is, in general, aliased into periods longer
than 20 days (Andersen, 1995a).
Though the orbit configuration for the TOPEX-class satellites was designed to
optimally recover ocean tides (Parke et al., 1987), the two largest semidiurnal
constituents, M2 and S2, are aliased into 62-day and 59-day period respectively. In
addition, at least ~3 years are required to separate M2 from S2. This serves as a critical
period to the design of the data time span for interleaved mission to ensure the
separability of M2 and S2 tides. Moreover, ~9 years (3355 days) of time series are
required to separate the K2/P1 and K1/Ssa pairs. Note that K1 is also aliased to the period
very close to semiannual signal (Ssa) (i.e. 173-day) that also causes the separation
problem between these two tides.
The aliased and Rayleigh periods for GFO and sun-synchronous Envisat satellites are
larger than that of TOPEX. That means more data are necessary to obtain reliable tidal
estimates at the same location. For GFO, M2 tide is aliased to the period close to the
period of annual signal (Sa) (i.e. 317-day). P1 tide is aliased to the period of ~12 years
(4464-day), however, the entire data time span of GFO lasted for 8 years only. Therefore,
it is not possible for GFO altimetry data to solve for P1 tide. It is also difficult to separate
M2/Sa tidal pair without 6 years of data time span. The tidal pairs for S2/K1, S2/Ssa, K1/ Ssa
also exhibit similar aliased period, as mentioned from the above.
Due to sun-synchronous orbit of Envisat satellite, S2 always has the same phase
during every repeat period. Hence, S2 tidal constituent appears as a constant height all the
time and it is thus impossible for Envisat to solve for. In addition, the M2/N2 tidal
constituent pair requires ~9 years (3169-day) of data time span for the separation. The
diurnal constituents K1 and P1 have the alias period of exactly 1 year (365-day). That
causes them inseparable from the annual signal (Sa) and from each other.
The aforementioned facts detail the requirement of minimum time span for a full
periodic tidal signal recovery and the difficulty in the separation among tidal constituents
using harmonic analysis at the same location. In contrast, the response formalism
(together with its orthogonal representation) does not have the above requirements
because the determined response weights are, in principle, independent of the Raleigh
criterion (Andersen, 1994). In other words, once the response weights have been
determined, a solution for any constituent can be inferred even for S2 constituent for sun-
synchronous satellite. However, the determined tides are not guaranteed to be of good
quality. Longer data time span is always required for better tidal estimates. Ray et al.
(2011) also argues that the assumption of smooth admittance across tidal band for the
response formalism is no longer valid when compound tides have significant amplitudes.
TOPEX altimetry data, though the presence of aliasing effect, serve as a backbone to
ocean tide determination because its tidal aliasing and Rayleigh periods are far better than
34
that for GFO and sun-synchronous Envisat satellites. The two-step empirical approach
based on along-track solution is not preferable, since the two major tidal constituents, M2
and S2, requires ~62 days to be solved and ~3 years to be separated from each other for
TOPEX mission. It also renders the ill-conditioned design matrix for GFO and Envisat if
a gridded tidal solution based on TOPEX alone is not initially determined to serve as a
constraint. In contrast, the one-step empirical approach acquires different altimetry
satellite data within a certain distance from the pre-defined grid center. This allows multi-
satellite data from different along-, adjacent-, and crossover- tracks to be utilized together
to achieve an empirical mitigation of the tidal aliasing effect.

3.3 Chapter Summary

We introduce the principle of satellite altimetry and illustrate the fundamental
difference in spatial and temporal sampling characteristics of satellite altimetry when
compared to tide gauges. Necessary instrumental, media, and geophysical corrections of
sea surface height (SSH) measurements are briefly described for the ocean tide modeling.
This provides a foundation for the linear observation setup for a weighted least-squares
solution in section 4.2.
We also introduce the tidal aliasing effect due to fundamental satellite orbit design
that the repeat period of altimetry satellites is far longer than that of the diurnal (daily)
and semi-diurnal (half-daily) tidal signals at a fixed location. This provides a better
understanding on the requirement of time span for the full tidal signal recovery and for
the separation among tidal constituents using harmonic analysis at the same location. The
response analysis, in contrast, does not have the above requirements since a solution for
any constituent can be inferred once the response weights have been determined, even the
S2 constituent from data of sun-synchronous satellite (e.g., Envisat). Nevertheless, more
data on temporal scale is always required for better tidal estimates.
TOPEX altimetry data is served as a backbone to ocean tide determination, no matter
which empirical approach to be used. To mitigate tidal aliasing effect, the one-step
empirical approach has its own merit in allowing multi-satellite data of different time
span from different along-, adjacent- and crossover- tracks to be utilized together.

35
 Chapter 4: Ocean Tides Modeling using Satellite Altimetry
In this chapter, we provide a review of ocean tides modeling approaches using
satellite altimetry. Based on the review, we illustrate the observation equation setup for
the ocean tide estimation on a regular grid, followed by proposing a spatio-temporal
combination approach that accounts for spatial and temporal (co-)variances
simultaneously for data weighting in a weighted least-squares solution process. Particular
emphasis is paid on the representation of spatio-temporal combination, and spatial and
temporal (co-)variance model specifications. Outlier detection criteria and robust
estimation technique are introduced for iterative re-weighting scheme.

4.1 A Review of Ocean Tides Modeling Approaches Using Satellite Altimetry

Satellite altimetry brought a new era in the study of ocean tides. It allows one to
generate a global ocean tide model from satellite altimetry based on equation (3.3) using
the tidal analysis methods described in Chapter 2. Three different approaches are applied
to generate ocean tide models. They are (1) hydrodynamic model; (2) assimilation model;
and (3) empirical model.
In the following subsection, some representative global ocean tide models, along with
their revised versions, are reviewed and classified based on their approaches and
methodologies.

4.1.1 Hydrodynamic modeling

The Laplace tidal equations (LTE), originally established by Laplace in 1775, and
discussed in section 2.3.1, is a dynamic tidal theory that describe the motion of ocean
water as a result of tidal forcing. These equations provide a solid foundation on the
dynamic principle of ocean tides. Due to the impossibility of obtaining an analytical
solution of the LTE, numerical methods have been the main driver to model ocean tides,
not to mention the non-linearity terms – such as bottom friction and advection –
incorporating into the dynamic equations for shallow water regions.

36
 Hydrodynamic models are derived by solving the LTE numerically, using bathymetry
data as input of depths, and ocean tidal constants observed by tide gauges around the
world as the boundary conditions (or data constraints). This modeling approach includes
simultaneously the solid-earth body tides, the ocean tidal loading, and the self-attraction
to solve for the ocean tidal height in the dynamic equations.
However, topographic drag and the bottom friction coefficient, describing the tidal
dissipation over the continental shelves and the shallow water regions, are of crucial
importance. In numerical modeling, empirical fine tuning of these two parameters in
specific regions are taken place to obtain a better solution, no matter realistic or not
(Matsumoto et al., 2000; Arbic et al., 2004). As a consequence, the tidal energy
dissipation is poorly simulated, in general.
Some models have treated this problem by using linear or quadratic parameterization
of bottom friction and by including the shallow areas in their domain of integration.
Some other models have treated this problem by assuming the ocean as frictionless, but
allowing energy to radiate through boundaries in the shallow water areas where energy is
dissipated.
One way to overcome this weakness is to increase the resolution. Another way is to
use the finite element method which improves the modeling of rapid changes in ocean
depth, the refinement of the grid in shallow waters, and the description of the
irregularities of the coastlines (Le Provost, 2001). However, a good bathymetry model is
required.
Two global hydrodynamic models have been derived in the past literature. The first
global numerical model was developed by Schwiderski (1980). It is constructed by
incorporating tide-gauge derived tidal constants and the best available bathymetry data
into the hydrodynamic interpolation scheme to solve for the LTE numerically. Despite
the dependence on the quality of the observations used and the existence of large errors in
this model, it had been the best available model for more than a decade until the era of
satellite altimetry. The spatial resolution of Schwiderski’s model is 1o  1o, except for
some semi-enclosed basins such as the Mediterranean Sea. It includes 8 major
constituents (M2, S2, N2, K2, P1, K1, O1, Q1) and 3 long-period tides (Ssa, Mm, Mf).
The other model, which was developed by Le Provost et al. (1994), is called FES94.1
model. This model is a modification of LTE that includes bottom friction parameterized
through a quadratic dependency on local tidal velocities, particularly for shallow water
regions. The equations were numerically solved by the finite element method (Le Provost
and Vincent, 1986). The FES94.1 model, with a resolution of 0.5o  0.5o, has a full
coverage of the world ocean including marginal seas and high latitudes, especially in
areas semi-covered by ice or under permanent ice shelves in the Weddel Sea and the Ross
Sea. This makes it as the default solution in most empirical models for the region beyond
 72o, due to limitation of the spatial coverage of TOPEX and GFO satellites. However,
this model is undefined in the Mediterranean Sea. It includes 8 major constituents (M2,
S2, N2, K2, P1, K1, O1, Q1) and 5 secondary constituents (Mu2, Nu2, L2, T2, 2N2), linearly
deduced by the admittance function along with nodal modulations and equilibrium long-
period tides.

37
4.1.2 Assimilation modeling
The assimilation modeling approach, which is essentially the hydrodynamic
approach, was proposed to integrate altimetry data and/or tide gauge data into the
hydrodynamic model for overcoming the weakness in Deep Ocean without data. Three
current global models are reviewed as an illustration in this modeling approach.

The NAO.99b is a global ocean tide model developed by Matsumoto et al. (2000). It
was based on the hydrodynamic tidal equations derived by Schwiderski (1980) with the
inclusion of advection terms and ocean tidal loading computation. It assimilated ~ 5 years
of TOPEX altimeter data into the hydrodynamic model. The sea surface height (SSH)
anomalies were analyzed through the response method with orthotide formulation. The
free-core nutation (FCN) resonance (Wahr, 1981 and Wahr and Sasao, 1981) 7 and the
radiational potential effect (Cartwright and Ray, 1994) 8 were taken into account in the
analysis. The tide model was provided on a 0.5o  0.5o grid.
The TPXO model, with the ongoing updated versions TPXO6.2, TPXO7.1, and the
latest TPXO7.2, was developed by Egbert et al. (1994) using the inverse scheme OTIS
(Oregon State University Tidal Inversion Software) to assimilate tide gauge and altimetry
observation data into the hydrodynamic equations by the representer approach (Egbert et
al., 1994, Egbert and Erofeeva, 2002). The eight major semidiurnal and diurnal tides are
provided together with two long-period tides (Mf, Mm) in the form of tidal harmonic
constants on a 0.25o  0.25o full global grid.
The FES model, with the ongoing update versions FES98/99, FES2002, and the latest
FES2004, is a series of models representing an improvement over its predecessor FES98
and FES94.1 (Lefèvre, 2000), in which only tide gauge data are assimilated into
hydrodynamic equations. In FES99, 670 tide gauges and 687 TOPEX altimetric
crossover data sets were assimilated by a revised representer method, similar to TPXO
models, to improve the accuracy for the FES98 model (King et al., 2005). For both
models, FES99 and FES2002, the solutions are provided on a 0.25o  0.25o global grid.
The latest model, FES2004, which assimilated ERS-1/-2 altimetry data in addition to
TOPEX, is provided on a 0.125o  0.125o global grid (Lyard et al., 2006).
It should be noted that both hydrodynamic and assimilation modeling approach
assimilate tide gauge data as the boundary conditions into hydrodynamic models to
compute ocean tides. Thus the evaluation using ground-truth tide gauge records is not
completely independent, for instance in case of FES2004 or the TPXO models. In
addition, the semi-empirical models (i.e., the GOT00.2/4.7, EOT08a/10a/11a, and
DTU10 model), which used these models as background model for higher spatial
resolution, would render the dependence on the ground-truth tide gauge records. As a

7
Free core nutation (FCN) resonance is a normal mode of the Earth, consisting of a relative rotation
between the fuild core and mantle together with associated deformation described by frequency-dependent
love numbers of degree 2.
8
Solar radiation is a factor that affects S2 and T2 tides indirectly. Therefore, this radiational effect in form
of radiational potential has to be included within the gravitational potential of S2 tide. This model
correction is done by rescaling the S2 amplitude with 0.97 and phase lagged by 5.9o.
38
result, it makes the external accuracy comparison with tide gauges not meaningful since
the models have already utilized the tide-gauge-derived harmonic constants.

4.1.3 Empirical modeling

The empirical modeling approach neither requires the knowledge of the bathymetry
or coastal geometry, the bottom friction coefficient and the topographic drag, and the
advection terms governing the tidal dissipation, nor the numerical scheme to solve for
hydrodynamic equations. It requires sufficient altimetry data available only for the ocean
tide determination. When both the hydrodynamic and the empirical modeling approaches
are inter-compared, the empirical modeling approach offers a simple, effective, and
practical method for tidal analysis and prediction.
Though it offers less physical insight and understanding when compared to the
hydrodynamic modeling approach, one can get into a more accurate representation of
ocean tides to quantify tidal dissipation and other related quantities. The hydrodynamic
and/or assimilation modeling approach also suffers from computational resource
requirement, due to its numerical nature, and the fact that inaccuracy arises from
inadequate bathymetry data and unknown friction and viscosity parameters (Ray et al.,
1996; Matsumoto et al., 2000).
As a result, many empirical ocean tide models are derived from satellite altimetry
data only, owing to the high precision altimeter measurements and better orbit
determination techniques since 1990s. Note that the altimeter measures geocentric tidal
height which includes ocean tide and ocean tidal loading, as well as solid Earth tide.
After correcting for solid Earth tide effect, one can solve for both ocean and ocean tidal
loading (i) in an iterative fashion (Ray, 1999), (ii) using the 7% rule of ocean tide for the
ocean tidal loading correction (Andersen, 1999), or (iii) the altimetry measurements are
previously corrected from forward loading models. Ocean tides are then determined via
harmonic analysis, response analysis (Munk and Cartwright, 1966) extended with
orthotide formalism (Groves and Reynolds, 1975; Cartwright and Ray, 1990b) or
Proudman functions (Sanchez and Pavlis, 1995).
There are two ways to generate the empirical models: (1) indirect analysis of
altimetry data from the background ocean tide model (called “semi-empirical method”)
and (2) direct analysis of altimetry data (called “purely empirical method”). In the first
method, the SSH is preliminarily corrected with an a-priori ocean tide model, followed
by using the incremental SSH anomaly for the incremental tidal solution, which is later
added back to the a-priori model to get the new full model. In the second method, a full
tide solution is derived by using the SSH anomaly derived from satellite altimetry
directly.

39
 Empirical models, no matter purely empirical or semi-empirical, are derived from two
methods. The first method is to conduct the along-track tidal analysis from TOPEX
altimeter data, followed by spatial interpolation onto a regular grid (Andersen, 1994,
1995; Desai and Wahr, 1995). Note that homogeneous weighting was applied to the
corresponding altimeter data. This amounts to the following two-step method; when GFO
and Envisat data are included, an iterative solution step from TOPEX data has to be
conducted. In other words, a predetermined TOPEX-alone regular gridded tidal analysis
solution is served as a background ocean tide model for the incremental tidal analysis of
GFO and Envisat altimeter data. This is because the tidal aliasing effect of GFO and
Envisat satellites are much worse than that of TOPEX, as discussed in section 3.2.
The second method is to acquire multi-satellite altimeter data at a certain distance
from a regular grid. This is followed by weighting the data spatially based on Gaussian
distance decay in the least-squares solution process (Eanes and Bettadpur, 1995; Smith,
1999; Savcenko and Bosch, 2008; Bosch et al., 2009). This amounts to a one-step
method. No temporal weighting has been applied. The rationale is the mitigation of the
tidal aliasing effect empirically by utilizing more along-tracks and cross-tracks of multi-
satellite altimeters with distinct spatial and temporal coverage. Note that most current
ocean tide models conduct the incremental tidal analysis based on FES2004 background
model.
The first altimetry-derived model was given by Cartwright and Ray (1990a, b, 1991).
This model was obtained through response analysis using the orthotide formulation,
based on 2.5 years of Geosat altimetry data. Since the launch of TOPEX in 1992, more
than 20 global tide models have been developed from the altimetry data. Considering the
limit of content, only five representative empirical models will be described as an
illustration of different modeling approaches in this classification.
The DW95 model is a purely empirical ocean tide model, developed by Desai and
Wahr (1995, 1996). The most updated version 7.0 was estimated from the observations
performed for the repeat cycles between 10 and 229 of the TOPEX altimeter mission.
The orthotide response formulation is utilized to solve for the diurnal and semidiurnal
ocean tides through admittance that linearly interpolated in spectral domain across
narrow bandwidths around each of the monthly (Mm), fortnightly (Mf), and termensual
(Mt) tidal components. This is a purely empirical model without reference to any a-priori
tide model or any direct or indirect information from the dynamics of the tides. The tidal
solution is estimated in grid size of 2.834o in longitude by 1o in latitude followed by
spatial interpolation onto 1o×1o grids within the limit of the TOPEX spatial coverage of
 66o. Beyond  66o, this model is extended with the Schwiderski’s model.
The CSR4.0 model is the revision of the older version CSR3.0 (Eanes and Bettadpur,
1995). The CSR ocean tide model series was developed by Eanes and others by using
TOPEX altimetry data and the orthotide formulation. The CSR4.0 model is obtained
using response analysis of about 6.4 years of the TOPEX altimetry data that solves for
incremental sea surface height (SSH) anomaly based on CSR3.0. Note that FES94.1
serves as a background model for CSR3.0. The corrections were produced in 3o  3o grids,
followed by spatial interpolation onto 0.5o  0.5o gridded resolution based on a

40
convolution with the 2-D Gaussian function. This model includes 8 major tidal
constituents.
The GOT model, with the ongoing update versions GOT99.2b, GOT00.2, and the
latest GOT4.7, was initially developed at NASA’s GSFC, and is known as SR94
(Schrama and Ray, 1994) and SR95.0/.1. Altogether, 232 cycles of TOPEX altimetry
data were used to derive the solution for 8 major semidiurnal and diurnal tides. The tides
were computed based on FES94.1 background model. The FES94.1 default model was
directly utilized outside the latitude limit of the TOPEX data (  66o). With better orbit
information available, all cycles of 10-day TOPEX data, complemented by 81 35-day
cycles of ERS-1/2 data in shallow seas and in polar seas, are used to solve for the
geocentric tide. This is followed by determining the ocean tides and ocean tidal loading
in an iterative manner. Also note that the a-priori models used in the GOT model include
not only FES94.1 model, but also some local and regional hydrodynamic models (Ray,
1999). The latest model, GOT4.7, is supplemented with more altimetry data and is
extended to include the S1 and M4 constituents, in addition to 8 major constituents. Minor
tides, which are inferred by the admittance function from the major tides, are also
available for computation with the software package. The tidal solutions are given on a
0.5o  0.5o grid.
The DTU10 model is a current global ocean tides model developed by the Technical
University of Denmark in 2010 (Cheng and Andersen, 2011). It was generated by using
the response method for an incremental analysis of multi-missions altimeter data based
on FES2004 background model. Phase A and four years of phase B data from multi-
mission altimetry measurements (TOPEX/POSEIDON and Jason-1/2), all cycles of the
altimeter datasets from ERS-2, Geosat Follow On (GFO), and Envisat were used at
latitude coverage of ±66o. The altimeter data are generated from the Radar Altimeter
Database System (RADS) database. Outside the coverage of these altimetry data, the
model was relaxed to the FES2004 default model. It includes 8 major semidiurnal and
diurnal tides, with the S1 and M4 constituents taken from the GOT4.7 model. The tidal
solutions are given on a 0.125o  0.125o grid. The model grid also extended onto the land
for proper accuracy assessment when compared against coastal tide gauges. The missing
grids were covered based on the DTU10BAT bathymetry land mask.
The EOT model, with ongoing update versions EOT08a, EOT10a, and the latest
EOT11a, was developed at the Deutsches Geodätisches Forschungsinstitut (DGFI).
EOT08a used entire cycles of TOPEX/Jason-1, GFO, ERS-1/-2/Envisat altimeter data
with refined orbit, new data pre-processing techniques, and improved geophysical
corrections, to derive the solution for 8 major tides, together with the 2N2 and M4
constituents. The tides were computed based on FES2004 background model. Different
strategies for the transition to polar oceans were separately applied to each tidal
constituent, as described in detail by Savcenko and Bosch (2008). The successor models,
EOT10a and EOT11a, utilized more SSH datasets for an incremental tidal analysis from
Jason-1/-2 and Envisat and excluded ERS-1 and GFO datasets, based on EOT08a as
background model. The data were upgraded by means of a newer orbit and improved
geophysical corrections for the incremental tidal solution based on the same background
model, followed by determining the ocean tides and ocean tidal loading in an iterative
41
manner (Savcenko and Bosch, 2010, 2011). The S1, Mm, and Mf constituents were
included in the solution, using harmonic tidal analysis, and the transition zone was
refined at latitudes higher than ±62o. The tidal solutions are given on a 0.125o  0.125o
grid.

4.2 A Novel Spatio-Temporal Combination Approach

In this dissertation, the one-step empirical method, which acquires multi-satellite
altimeter data at a search area from a predefined regular grid center location, is employed.
Each grid was set up at a constant 0.25o interval in both latitude and longitude direction.
A search area for acquiring multi-satellite altimeter data is defined as a 3×3 data window,
corresponding to 0.75o× 0.75o for the tidal analysis at the grid center ( c , c ) (Figure
4.1).
Since altimetric satellites repeat sampling the sea surface height (SSH) at the same
location ( ,  ) every 10-day for TOPEX, 17-days for GFO, or 35-day for Envisat,
respectively, each location along the satellite altimetry ground track contains the time
series of SSH anomaly. The tidal solution at the grid center, ( c , c ) , is then estimated
using these SSH anomaly time series at different locations via weighted least-squares in
one-step. Note that the time series at different locations is not first reduced a-priori to the
grid center, but is reduced and solved simultaneously at once in a weighted least-squares
process.

Figure 4.1 Schematic diagram of the search area to acquire TOPEX (Blue), GFO (Yellow), and Envisat
(Red) altimeter data for different satellite ground tracks with an angular distance
d  ((  c ) cos ) 2  (   c ) 2 from each predefined regular grid center location (c , c ) .

42
 Moreover, different satellite altimeters have different data time span and ground
tracks. For example, TOPEX altimeter data were collected between Aug 1992 and Dec
2005, whereas that of GFO was collected between Jan 2000 and Nov 2007 (Table 3.1 for
details). The SSH measurements at different time epochs are sampled independently.
Similar situation applies to that in space because the sampling locations from different
satellite passes for each altimeter satellite are pre-defined due to satellite orbit design.
Hence, no covariance information is available among different altimeter datasets at
different time epochs and at different locations. However, empirical temporal and spatial
correlation relationships do exist and can be constructed as will be discussed in section
4.2.2, but it is not considered in this dissertation as a matter of time efficiency and
simplicity.

4.2.1 Observation Equation Setup

Figure 4.1 defines the basic setup for the observation equation. The time series of
SSH anomaly observations of each location in space, ( , , t ) , for different satellite
altimetry ground tracks are used to estimate the tidal solution at the grid center location
( c , c ) . An appropriate (co-)variance matrix for data weighting both in space and in
time simultaneously has to be constructed, instead of just equally weighted or weighting
in space that have been applied in the previous literature for one-step empirical method.
Note that the weight matrix usually comes from the inversion of the covariance matrix in
geodetic science, but exceptions are not uncommon in geographical (Fotheringham et al.,
2002), geodetic, and geophysical applications (Junkins et al., 1973; Guo et al., 2006;
Sabaka et al., 2010).
The orthotide formulation, an orthogonal representation of the response method for
the eight major tides proposed by Groves and Reynolds (1975), is adopted to estimate the
unknown coefficients of orthotide functions. A modification to the orthotide formulation
is made to include satellite pass biases, trend, seasonal signals (i.e., annual Sa , semi-
annual S sa ), long period tidal components (i.e., monthly M m , fortnightly M f ), and short
period tidal components (i.e., M 4 and S1 ). Therefore, the observed instantaneous SSH
anomaly at time t, hssha obs
( ,  ,t ) , without ocean tidal correction is expressed in an
extended orthotide formulation as:
obs
hssha ( ,  , t )  a i ( c , c )  b( c , c )(t  t 0 )
2 2
  [U mj ( c , c ) Pjm (t )  V jm ( c , c )Q mj (t )]
m 1 j  0
(4.1)
6
  [ f k C k ( c , c ) cos( k (t )   k   k )
k 1

 f k S k ( c , c ) sin( k (t )   k   k )]  e( ,  ,  c , c , t )

43
where Pjm (t ) and Q mj (t ) are the orthotide functions calculated from the tide generating
potential, t0 is the reference time epoch,  k (t ) are the Doodson arguments at Greenwich
for each constituent,  k are the additional phase corrections, f k and k are nodal factors
and angles, respectively, applied to the lunar constituents (i.e., M 4 in this case), and
e( ,  ,  c , c , t ) includes the random noise, and the errors due to geophysical and media
corrections, the altimeter temporal measurement noise at location ( ,  ) due to the
general ocean circulation signal at time epoch t, and the error due to spatial variation of
the general ocean circulation around the pre-defined grid center location ( c , c ) since
tides are not estimated at the observed locations but at ( c , c ) .
At least 26 unknown parameters are to be solved, depending on the number of
satellite altimeter passes within the search area. The unknown parameters to be estimated
in equation (4.1) are listed as follows:
(1) ai ( c , c ) is the local bias of satellite altimeter pass i;
(2) b( c , c ) is the local sea-level trend;
(3) U mj ( c , c ) and V jm ( c , c ) are the unknown location-dependent coefficients of
orthotide functions, which can later be converted back to amplitude and phase of major
tidal constituents through equation (2.45) and (2.49);
(4) Ck ( c , c ) and S k ( c , c ) are the unknown in-phase and quadrature amplitudes
for the four long-period tides (i.e. Sa , S sa , M m , M f ) and other short-period tides (i.e.,
M 4 , S1 ).
In this study, both spatial and temporal (co-)variances will be taken into account
simultaneously in a balanced manner for data weighting of the tidal solution introduced
in the following sub-section. They have to be defined according to the underlining spatial
variation of the general ocean circulation around grid center, ( c , c ) , and the overall
temporal measurement noise (including the altimetry measurement noise due to the
general ocean circulation and an overall error budget including other errors due to
geophysical and media corrections) at a particular time t, respectively.

4.2.2 Weighted Least-Squares Solution, Representation of Spatio-Temporal

Combination, and Covariance Model Specification

Weighted Least-Squares Solution

Though the sea surface height (SSH) anomaly measurements still contain the
remaining errors due to the general ocean circulation signal, the error in geophysical and
media corrections, and the random noise (as discussed in section 3.1.2), the response tidal
analysis formulation is assumed to be linear; therefore, let us consider the Gauss-Markov
model (Koch, 1999) which is defined as follows:
44
 y  A  e , e ~ (0,    02 Qe ) with Q y  Qe : P 1 (4.2)
where y is the n × 1 observation vector, collecting the observed instantaneous SSH
obs
anomalies without ocean tides correction at each epoch time t, hssha ( ,  , t ) , affected by
the n × 1 random error vector e, A is the n × m coefficient matrix with full column rank,
 is the (unknown) m × 1 parameter vector to be estimated, such as the satellite pass
biases, sea-level trends and coefficients of orthotide functions, P is the symmetric and
positive-definite n × n weight matrix (still to be specified), and  02 is the variance
component of unit weight, with cofactor matrix Q y  Qe : P 1 identical to the inverse
weight matrix.  is the covariance matrix that is proportional to Q y or Qe by  02 .
Based on the principle of e T Pe  min , the weighted LEast-Squares Solution (LESS)
of the unknown parameter vector  is given by
ˆ  N 1 c , with N , c : AT P A, y  (4.3)
with the corresponding residual vector and its cofactor matrix given as
e~  y  Aˆ  Ry  (Q~e P) y  Re (4.4)
Qe~  Q y  AN 1 AT  RQe R T (4.5)
where R  I n  AN A P   Q~e P . The matrix R, sometimes called the redundancy
1 T

matrix, contains useful information in the sensitivity analysis for the effects of outliers
(Huber, 1981; Ding and Coleman, 1996), with the property that
n
r : tr ( R )   ri  n  m , for 0  ri  1 if P is diagonal, (4.6)
i 1

represents the total redundancy or the degree of freedom (r) of the model, where ri is the
i-th diagonal element of R, called the redundancy number, that indicates the good control
or high reliability when redundancy numbers close to 1 (Schaffrin, 1997). The estimated
variance component is given by:
~
e T P~
e
ˆ 02  (4.7)
r
which indicates the goodness-of-fit between the model and the observations.

Representation of Spatio-Temporal Combination

Contemporary empirical ocean tide models are estimated by equally weighted
solution or spatially weighted solution based on spatial (co-)variances. However, because
of the availability of the space-time data, it is justifiable that both spatial and temporal
(co-)variances has to be taken into account simultaneously, no matter for spatio-temporal
process prediction or for data weighting in the gridding process. The latter is the main
theme of this dissertation study. Hereafter, we refer the data weighting approach to as
spatio-temporal combination. A brief literature review about the development is
described, and the spatio-temporal combination approach is then proposed.
Gridding is a process to convert spatially (and temporally) scattered individual data
points with distinct random errors into values at regular grids. Simple or Ordinary
45
Kriging (equivalent to least-squares collocation (LSC) or the representer method), spline
minimum curvature, nearest neighbor, polynomial regression, expression using radial
basis function such as multi-quadrics, and triangulation are common methods, principles,
and algorithms that can deal with spatial and/or temporal data for the aforementioned
purpose.
Both the spatial and temporal processes were separately developed at the beginning of
1970s, since the data are obtained either in space or in time (Kyriakidis and Journel,
1999). The spatial processes, particularly in the field of geophysics and geodesy (Jordan,
1972; Moritz, 1978), were modeled based on covariance functions for least-squares
collocation, whereas the temporal processes were usually modeled using time series
analysis technique, such as autoregressive and moving-average processes (Box and
Jenkins, 1976).
Thanks to the advance of satellite geodesy, a large amount of space-time datasets are
available. Considerable attention has been paid to analyzing spatio-temporal phenomena
as a whole recently (De Cesare et al., 2002; Cheng, 2004). The fundamental concept of
analyzing spatio-temporal phenomena is to express a spatio-temporal covariance function
(with spatial lag s and temporal lag t ) in a product (Rodriguez-Iturbe and Meija,
1974) or an addition (Rouhani and Hall, 1989) of pure spatial and purely temporal
covariance functions; and hence, the pure spatial and temporal covariance matrices. The
expression of spatio-temporal covariance function in product model has been recently
applied to the mapping of the spatio-temporal sea surface height (SSH) anomaly residuals
(Le Traon et al., 1998).
However, it is impossible to tell which expression – the addition or the product form
– is a better representation when the exact process is unknown. To accommodate both
expressions, a product-sum modeling concept was proposed (De Cesare et al., 2001,
2002). Note that both spatial and temporal covariance models are being used for
prediction. However, in this dissertation, the concept is merely utilized to develop the
spatio-temporal combination approach for data weighting. Note also that the temporal
(co-)variances in this dissertation may not be entirely stationary because the measurement
noise level time epochs t and t  t may be entirely different depending on the sea state
governed by the general ocean circulation signal. This will be described later in this
section.
In this dissertation, the spatio-temporal combination approach is to define the
covariance matrix,  , to weight spatio-temporally obtained altimetry data in the
weighted least-squares process. The covariance matrix, based on the conventional
covariance matrix setup in geodetic science, can be expressed as (De Cesare et al., 2001,
2002):
   S  T   ST
(4.8)
  S2 QS   T2 QT   ST
2
QST
where QS , QT , and QST : QS  QT   [QS ij QT ij ] are the n  n observation cofactor
matrices for the spatial (co-)variance, the temporal (co-)variance, and the Hadamard
product for the multiplication of the spatial and the temporal covariance function in

46
matrix form, respectively. The Hadamard product results in a symmetric and positive-
definite matrix.  S2 ,  T2 , and  ST
2
are the scaling factors to be determined, called variance
components, for space only, time only, and their product, respectively. The modeling
solution employing this covariance matrix with specific estimates ˆ S2 , ˆ T2 , and ˆ ST2
for
data weighting is referred to as OSU12vce model.
Besides the conventional covariance matrix setup in geodetic science in equation
(4.8), an empirical data weighting scheme is not uncommon in the form of a conjectured
weighting function or a set of ad-hoc procedures whenever a better result can be obtained
and/or no precision information is provided for the description of data quality. This
scheme has been applied in geographical (Fotheringham et al., 2002), geodetic, and
geophysical applications (Junkins et al., 1973; Grafarend et al., 1980; Guo et al., 2006;
Sabaka et al., 2010). Based on such a conjecture, another covariance matrix,  , may be
written as:
1
   
   S2 PS  T2 PT  ST2 PST    02 S PS  T PT  ST PST 
1
(4.9)
0 0 0 
where PS  Q S1 , PT  QT1 , and PST : QST
1
are the inverse of n  n cofactor matrices for
the spatial (co-)variance, the temporal (co-)variance, and the Hadamard product,
respectively while S , T , and  ST are known as “weight components”. The model
solution employing this covariance matrix with specific estimates ̂ S , ̂T , and ̂ ST
refers to OSU12sw model.
The cofactor matrices QS and QT will be specified below. The techniques of
determining the variance components for the equation (4.8) and the weight components
for the equation (4.9) will be illustrated in section 4.2.3.

Covariance Model Specification

To completely define the covariance matrix,  , for data weighting, both the temporal
and spatial (co-)variances have to be specified. The temporal (co-)variance can be
specified by temporal noise and assumption of serial correlation structure. We describe
more in detail on the spatial (co-)variance modeling because tides are prone to localized
effects that are quite different even when only a few kilometers apart, particularly near
coastal regions and over continental shelves, in order to improve the ocean tides along the
world’s coastal regions and over continental shelves when compared to existing models.
Note that the cofactor matrices for spatial (co-)variance, QS , and for the temporal (co-
)variance, QT , are assumed to be diagonal throughout the study, as mentioned, for time
efficiency and simplicity because matrix size increases with the number of observations
that poses difficulties for the matrix inverse and for matrix multiplications in the
weighted least-squares solution process, even for the block-diagonal matrices.
For temporal (co-)variances, we take into consideration the noise level of each 1-Hz
along-track measurement for each altimeter plus the total error budget from different
error sources. Significant wave height (SWH) is a quantity that describes the wind-driven
47
wave heights, calculated from the standard deviation of the highest 30% of the measured
sea-surface wave displacement (Stewart, 2007). The value of SWH, which is one of the
measurements of radar altimeter, depends on the sea state of a fixed location at a
particular time epoch t. In other words, if the sea state of the location is not stable at a
time epoch t, the standard deviation of the measured sea-surface wave displacement is
high, and hence, the SWH value.
Figure 4.2 depicts approximately linear relationship between SWH and altimeter
noise (Fu et al., 1994). This relationship can be used to recover the altimeter noise level
approximately from the SWH measurement at a different time epoch t. The sum of the
variance of the altimeter noise, the variance of random orbit error of each satellite, and
the variance of other random errors form the variance of the total noise,  t2 , at time t for
the temporal variance. Since each measurement for each time epoch is independent, the
off-diagonal element (i.e. covariance) is assumed to be zero in this dissertation for
simplicity. Hence, the diagonal of the cofactor matrix for the temporal variance, QT ii , is
formed as  t2 , so that heterogeneous temporal variances are considered and incorporated,
instead of homogenous one. An overview of the major error components in the altimeter
measurement corrections and SSH measurements nowadays are listed in Table 4.1.

Figure 4.2 The relationship between altimeter noise and Significant Wave Height (SWH) (Adopted from Fu
et al., 1994)

Table 4.1 Standard deviation (±cm) of sea surface height (SSH) measurements of TOPEX/POSEIDON,
GFO and Envisat

Major error sources TOPEX/Jason-1/-2 GFO Envisat

Orbita 2.0/1.5/1.5 5.0 2.0
Ionosphereb 0.5 0.5 0.5
Troposphereb 1.2 1.2 1.2
Other errors 2.0 2.0 2.0
Altimeter noisec 1.6-3.0 2.0-5.0 1.8-2.3
a. Reference: TOPEX/Jason-1/-2 (Bosch, 2004; Menard et al., 2003; Desai and Haines, 2010), GFO
(Lemoine et al., 2006), and Envisat (Rudenko et al., 2012), respectively.

48
b. Reference to Menard et al., 2003.
c. Altimeter noise range for Significant Wave Height (SWH) between 2 m and 5 m.

For spatial (co-)variances, we take into account the rapid distance decay property near
coasts in Shallow Ocean. In other words, for the altimeter measurements very close to the
pre-defined grid-center estimate, a very small variance should be assigned for the
gridded tidal estimates, in particular very close to the coast to preserve the locality of
tides. This is because tidal dynamics, changes significantly when the tidal wave enters the
continental shelf and shallow water regions from the deep oceans as manifested from the
tidal wavelengths. Our method significantly modifies the method described by Andersen
(1999).
Andersen (1999) proposed an isotropic Markov covariance function (i.e. an
exponential model) that incorporates the tidal wavelength property of diurnal and semi-
diurnal waves for the least-squares collocation of the obtained along-track tidal solutions
when entering shallow water regions, since two-step empirical method is employed. The
presence of a highly varying water depth also changes the tidal wave property near the
coasts. The covariance function below is utilized for the least-squares collocation
estimates of each tide reduced from the along-track tidal solution to the grid
center ( c , c ) :
d
d ( )
C (d )  C0 (1  )e  (  m )
 ( m ) (4.10)
where C0 is the error variance,  ( m ) is a correlation length (also called characteristic
length) which is a function of the tidal wavelength, d  ((    c ) cos  ) 2  (   c ) 2 is
the distance between the observations at (, ) and the grid center ( c , c ) , as defined in
Figure 4.1, and  m is the tidal wavelength for each tidal species m which is related to the
bathymetry depth H b (Pugh, 1987) for the correlation length which is given by:
 m  ( gH b ) 1 / 2 Tm (4.11)
where g is the acceleration due to gravity and Tm is the period of each tidal species m.
The bathymetry depths versus tidal wavelengths for diurnal, semi-diurnal, and quarter-
diurnal constituents are tabulated in Table 4.2. The value of correlation length governing
the decay is empirically defined based on the range of bathymetry depth in shallow water
(Andersen, 1999). Similar but different empirical values of correlation length have
currently been implemented in the DTU10 model (Cheng and Andersen, 2011).

Table 4.2 Tidal wavelength versus depth for different tidal species.

Wavelength (km)
Depth (m)
Diurnal Semi-diurnal Quarter-diurnal
5000 19811 9905 4953
1000 8860 4430 2215
500 6265 3132 1566
100 2802 1401 700
49
 10 886 443 221

In contrast, spatial (co-)variances are specified differently for data weighting in one-
step least-squares solution process in this dissertation instead of gridding along-track tidal
estimates by least-squares collocation in two-step. The diagonal of the cofactor matrix for
the spatial variance, QS ii is formed as,
s 2,
QS ii  (4.12)
w(d , Λ) cos
where s2, is the location-dependent variance of general ocean circulation signal as
shown in Figure 3.3, cos  is the function taking into account the latitude convergence of
the poles, and w(d , Λ ) is the function that takes into consideration of the tidal
wavelength changes in shallow water regions via Gaussian distance decay property for
the diurnal, semi-diurnal, and quarter-diurnal tidal species included in the observations.
Hence, the longer the distance between the observation locations and the grid center, the
larger the variance. Because the tidal species (i.e., diurnal, semi-diurnal, and quarter-
diurnal) can be independently separable through harmonic analysis or other means, it can
be expressed in the product form as,
2
 
3
d     m   m 
2  d 

w(d , Λ)   1    m   e (4.13)
   m  
m 1
 
where  m is set to 1, 2 , and 2 to account for half the decrease in tidal wavelength from
diurnal to semi-diurnal, and from semi-diurnal to quarter-diurnal species (Table 4.2), and
 is an empirical value that scale all  m , with m =1 (diurnal), m =2 (semi-diurnal), and
m =3 (quarter-diurnal), with Λ   1  2  3  is the vector for the tidal wavelength of
the three tidal species. Note that the bathymetry depth in this dissertation is extracted
from ETOPO1 global relief model, which is a 1 arc-minute model of Earth’s surface that
integrates land topography and ocean bathymetry (Amante and Eakins, 2009).
The characteristics of spatial (co-)variances used in this study are illustrated in Figure
4.3. It demonstrates rapid inflation of variances when entering into the shallow ocean.

50
Figure 4.3 Characteristics of spatial (co-)variances for the equation (4.5) under different depths for   40
and s2,  100 cm 2

In fact, possible empirical covariances, for the block-diagonal cofactor matrix QT or

for QS , can be incorporated when both temporal and spatial correlations are substantial.
For the temporal cofactor matrix, QT , the noise in the data time series may be assumed to
have a serial correlation structure. The estimated covariances from the autoregressive
coefficient can then be assigned to the off-diagonals of the cofactor matrix for the
weighted least-squares solution (Cochrane and Orcutt, 1949). It has been applied to
observed VLBI rates (Iz and Chen, 1999) and tide gauge data analysis (Barbosa et al.,
2008). However, it failed to get a better solution when compared to ordinary least-squares
(OLS), even for measurements with higher sampling rate when compared to that of
satellite altimetry, and for a magnitude of autoregressive coefficient estimate up to 0.6 in
the above two contributions. This can be successful only when the magnitude of
correlation coefficient is close to 1 (Iz, 2008).
An attempt for this method has been made. However, the magnitude of the
autoregressive coefficient is less than 0.08 in our altimetric tidal analysis in the Alaska
sea region with slightly inflated estimated standard deviation in the tidal estimates, when
compared to that of weighted LESS. The small magnitude of the autoregressive
coefficient is attributable to the long repeat periods of satellite altimetry (9.9156-day for
TOPEX/Jason-1/Jason-2; 17.0505-day for GFO; 35-day for Envisat).
51
 For spatial cofactor matrix, QS , the locations of the altimeter ground tracks were
fixed due to the fundamental orbit design. Arrangement in block-diagonal matrices is
possible provided that a pre-defined spatial covariance function is an appropriate
representation. An effort should be paid to data re-structuring. Note that the block
diagonal or full covariance matrices would increase the computation time as compared to
diagonal matrix due to a series of multiplications of several matrices along with their
matrices’ inversion. An approximation algorithm based on Monte Carlo’s simulation can
lessen the amount of computation time (Kusche, 2003), particularly for the variance
component estimation technique which will be described in the below section. But how
much the gain in time efficiency is still a matter in practice.
Note also that one requires the estimation of ~ 500,000 gridded tidal solutions for the
global ocean tide model at a 0.25o×0.25o resolution, with number of observations per
gridded tidal solution vary from 80 (e.g., for coastal regions where an estuary system is
present) to 35000 (e.g., at latitude higher than 60o in deep ocean) in this dissertation. This
is a computationally demanding task. A feasibility study of incorporating block-diagonal
empirical covariance matrices has to be conducted for altimetry data before the operation.

4.2.3 Techniques for Spatio-Temporal Combination

Two techniques for the spatio-temporal combination approach are proposed. The first
technique is the application of variance component estimation (VCE) of type repro-
BIQUUE that correlates the spatio-temporal multisatellite altimeter data in a balanced
manner with the covariance matrix specified in equation (4.8) for data weighting. The
second technique is a set of empirical procedures to specify the conjectured weight
components in equation (4.9).
Since the advancement of satellite geodesy together with ground-based observations,
it brings about heterogeneous data sets with different precision information no matter
whether from satellites or other sensors. These data sets have to be properly correlated
when they are combined in order to obtain the best information on the parameters from
weighted least-squares solution.

Variance Component Estimation (VCE) technique

In geodetic science, variance component estimation originated from Helmert (1907,
1924); it is a statistical technique that uses least-squares residuals to estimate
heterogeneous variance components, usually by iteratively reweighting the heterogeneous
datasets. A review of different types of VCE can be found, for example, in Schaffrin
(1983), Fotopoulos (2003), and van Loon (2008). It has been applied gyrocompass
azimuth observation variance analysis (Kleusberg and Grafarend, 1981). Contemporary
applications include GPS position time series, network adjustment and analysis
(Schaffrin and Iz, 2001; Yang et al., 2005; Amiri-Simkooei, 2009), regional
GPS/leveling/geoid network unification (Kotsakis and Sideris, 1999; Fotopoulos, 2003),
and satellite gravity data analysis (van Loon, 2008). It is anticipated that the application
of VCE theory is continuously growing given multiple data sources from different

52
satellites, airborne measurements, and ground-based stations, along with better
computational resources.
The fundamental concept of a Variance Component Model (VCM) is to define the
stochastic model, (i.e. the variance-covariance matrix)  in equation (4.8), with the
following covariance structure
  12 C11  12 C12   1 C1 
 
  12 C12T  22 C 22   2 C 2 
 . (4.14)
     
 
 1 C1  2 C 2    C 
T T 2

The partition of the matrices, Cij, is based on different distinct groups of observation
datasets. Then, the (co-)variance components,  12 ,  12 ,...,  2 , which are collected in the
(  1) / 2 1 vector  , are to be determined. Hence, the general structure is given by
(see, e.g., Grafarend et al., 1980; Schaffrin, 1983),
  1 
    i2 Qii    ij Qij , for i  j , (4.15)
i 1 i 1 j  2

with
 0  0  0  0
 
    
 0  0  C  0
0  0  ij

 
Qii :   C ii   and Qij :       (4.16)
0  0  T 
   0  Cij  0  0 
    
 
 0  0  0  0
For most applications in geodetic science, different groups of observation datasets are
uncorrelated with each other, such that the covariance matrix can be simplified as

    i2 Qii (4.17)
i 1

In this case, the Helmert estimates of the variance component, collected in the  1vector
   12  22   2  , are obtained from
T

Hˆ  q (4.18)
where
H ij  tr (Wˆ QiiWˆ Q jj ) and qi  e~ T (ˆ 1Qii ˆ 1 )~
e (4.19)
 
with Wˆ  ˆ 1 I n  A( AT ˆ 1 A) 1 AT ˆ 1 . Let 0 be the initially assigned covariance
matrix, then iteration is performed in equation (4.3), (4.4), (4.18), and (4.19) to estimate
ˆ which is used to determine ̂ until no elements in the vector ˆ of the estimated
variance components is changed in the k-th iteration (i.e., ˆ ( k )  ˆ ( k 1) ); this solution

53
represents the reproducing Best Invariant Quadratic Uniformly Unbiased Estimates
(repro-BIQUUE) according to Schaffrin (1983). Note that, in this case of uncorrelated
observation groups, the above VCM of the Helmert type also represents MInimum Norm
Quadratic Unbiased Estimates (MINQUE) (Grafarend et al., 1980; Koch, 1999). If the
different observation groups are correlated, it is neither of Helmert type nor MINQUE,
just the repro-BIQUUE. The number of iterations required will depend on the data
precision and the initial values.
Similarly to the conventional practice for iterative reweighting of the heterogeneous
datasets, we apply the above VCE procedures to the spatio-temporal combination
approach in equation (4.8). Note that Qii in equation (4.17) can be a full cofactor matrix
without the loss of generality. More explicitly, the equation (4.8) can be rewritten in
diagonal matrix form as,
   S2 QS   T2 QT   ST
2
QST
 s S21 0  0  sT21 0  0
   
2 0 s S2 2  0 2 0 sT2 2  0
S T
           
 2   2 
 0 0  s Sn  nn  0 0  sTn  nn
 s ST
2
0  0  (4.20)
1
 2 
2  0 s ST  0 
  ST 2
     
 2 
 0 0  s STn  nn
Helmert’s VCE technique, in the case, is to obtain ˆ  ˆ S2 ˆ T2  iteratively by
T
ˆ ST
2

equation (4.18) with equation (4.19) in the form of:

 tr (Wˆ Q S Wˆ Q S ) tr (Wˆ Q S Wˆ QT ) tr (Wˆ Q S Wˆ Q ST ) 
 
H   tr (Wˆ QT Wˆ Q S ) tr (Wˆ QT Wˆ QT ) tr (Wˆ QT Wˆ Q ST )  (4.21)
tr (Wˆ Q Wˆ Q ) tr (Wˆ Q Wˆ Q ) tr (Wˆ Q Wˆ Q ) 
 ST S ST T ST ST 

and
 e~ T (ˆ 1Q S ˆ 1 )e~ 
 
q   e~ T (ˆ 1QT ˆ 1 )e~  (4.22)
 ~ T ˆ 1 ˆ 1 ~ 
 e ( Q ST  )e 
 
where Wˆ  ˆ 1 I n  A( AT ˆ 1 A) 1 AT ˆ 1 . The initial covariance matrix, 0 , is assigned
by setting       1 . Then the iteration is performed in equation (4.3) and (4.4)
2
S
2
T
2
ST

for equation (4.21) and (4.22) iteratively until the estimated variance components in
vector ˆ  ˆ S2 ˆ T2 ˆ ST
2 T
 are no longer changed in the k-th iteration (i.e., ˆ (k )  ˆ (k 1) ).
The model solution employing this VCM is referred to as OSU12vce model. This
technique should provide a theoretical basis for optimal spatio-temporal data weighting.
54
Empirical procedures for conjectured weight component estimation
Besides the conventional covariance matrix setup in geodetic science in equation
(4.20) and the theoretically sound technique for iterative re-weighting, we describe a set
of empirical procedures that estimates the weight components in equation (4.9) and that
also generates a good result in this dissertation.
Because the set of observations is obtained in space-time, the cofactor matrices QS ,
QT , and QST fundamentally have the same dimension in this dissertation. Note that
contemporary empirical ocean tide models are achieved by equally weighted solution or
spatially weighted solution based on spatial (co-)variances. One can also conduct an
experiment on the trial of any cofactor matrix for data weighting. Our empirical
procedures are as follows:
Step 1: Using the equations (4.3), (4.4), and (4.7) within three separate models, one
can estimate the variance components,
e~ P e~ e~ P e~ e~ P ~
T T T
e
~S2  S S S , ~T2  T T T , and ~ST 2
 ST ST ST
r r r
independently and respectively, from the trial of n  n cofactor matrix for the spatial
(co-)variance, Q S  PS1 , the temporal (co-)variance, QT  PT1 , and the Hadamard
product, Q  P 1 , via a weighted least-squares process. Note that e~ , e~ , and ~
ST ST S T e ST

indicate the residual vectors that are predicted based on QS , QT , and QST , respectively.
The approximated weight components ˆ S  ~ S2 , ˆ T  ~T2 , and ˆ ST  ~ ST
2
in this
procedure are obviously different from inverse the variance components ˆ S2 , ˆ T2 , and
ˆ ST
2
estimated by the VCE technique above.
Step 2: The approximated weight components, ̂ S , ̂T , and ̂ ST are then used to
specify the covariance matrix,  , in equation (4.9). Then, this conjectured covariance
matrix is used for the data weighting in the weighted least-squares solution. The model
solution employing the empirical procedure refers to OSU12sw model.
To examine the effectiveness of both OSU12vce and OSU12sw model solution, a
comparative analysis is conducted to test against the equally weighted solution, the
weighted solution based on spatial (co-)variances, and the weighted solution based on
temporal (co-)variances with the same model specification, which will be shown in
Chapter 5.
In the next sub-section, we will describe the robust estimation technique for outlier
detection and down-weighting.

4.2.4 Outlier Detection and Robust Estimation

In both pre-processing and post-processing steps for any kind of observations, data
outliers – the observations inconsistent with the rest of the data set (Barnett and Lewis,
1994) – are ubiquitous and have a significant influence on the geodetic parameter
55
estimation (Krarup et al., 1980; Jørgensen et al., 1985; Kubik et al., 1985; Rangelova et
al., 2009). Even for a small number of undetected outliers (<0.2% of all data points), the
influence is substantial (Kern et al., 2005). They have to be detected and removed, or
controlled to ensure minimal distortion of the final parameter estimates.
The control of the influence of the outliers falls into three categories: (i) the empirical
procedures by setting thresholds such as the upper/lower bound with respect to the mean
(Barnett and Lewis, 1994) or median value (Hampel et al., 1986; Rousseeuw and Leroy,
1987), and the maximum allowable rates of change from one value to the next (e.g., Kern
et al., 2005); (ii) the conventional outlier detection test procedures based on statistical
testing (e.g., Baarda 1968; Pope, 1976) followed by their removal; and (iii) the robust
estimation techniques (e.g., Andrews, 1974; Krarup et al., 1980; Huber, 1964, 1981;
Hampel et al., 1986; Yang, 1999). The first category is usually utilized in a pre-
processing step for preliminary data screening for apparent and large blunders, whereas
the latter two are performed during or after the geodetic parameter estimation process in
an iterative fashion.
While the conventional outlier detection test procedures provide a significant insight
into both theoretical and practical considerations, the custom of detecting one susceptible
error at a time may render this procedure ineffective. It is because the resultant residuals
are sensitive to and potentially biased by outliers in other observations (Yang 1999; Guo
et al., 2010) due to the minimization of the sum of weighted squared residuals within a
Gauss-Markov (or any other) model. As a result, this causes good observations to be
flagged (Baselga, 2007), and potential loss of information when data points are
irregularly distributed over large areas with minimal data for the estimation process
(Rangelova et al., 2009). Some of these effects can be reduced, however, by following
the strategy of Koch (1983). On the other hand, they are also computationally expensive
when multiple outliers are suspected to exist for a large scale LESS problem (Kern et al.,
2005). In contrast, the robust estimation techniques provide less influence on the
estimates. This was shown by several outlier simulation studies for classical geodesy
problems, such as geodetic network analysis and GPS positioning (Jørgensen et al., 1985;
Hekimoglu and Erenoglu, 2007; Knight and Wang, 2009; Sisman, 2010). Given the huge
amount of data acquired both in space and in time since the advance of satellite geodesy,
the application of robust estimation technique in geodetic parameter estimation is slowly
gaining attention in the geodetic community (Krarup et al., 1980; Awange and Aduol,
1999). This technique is particularly helpful when a large-scale LESS problem is solved
and multiple outliers are presented, as long as a diagonal cofactor matrix can be assumed.

Outlier detection
Conventional practices include Baarda’s data snooping, Pope’s tau test, t-test, and F-
test. These tests use the normal, the τ, the t or the F- statistic to test against a critical or
threshold value for outlier detection and down-weighting certain observations in the
robust estimation. The above-mentioned methods should converge asymptotically by the
central limit theorem (e.g., Johnson and Wichern, 1992) when a large sample of
observations is involved and provided that observations are independent.

56
 Arranging the i-th residual in standardized or studentized form yields the test statistic
for data-snooping (Baarda, 1968), the τ (Pope, 1976), the t (student), and the F (Fisher)
test, with critical values of their corresponding distributions respectively as:
e~i
Wi   N  (0,1) (4.23)
 0 Qe~ ii 2

e~i 2
Wi  2
2
  12, (4.23a)
 0 Qe~ ii
e~i
i    (4.24)
ˆ 0 Q~e ii r,
2

e~i
ti  t  (4.25)
ˆ 0i Q~e ii r 1,
2
~
e 2
F : ti2  2 i  F1, f 1, (4.26)
ˆ 0i Qe~ ii
under the null hypothesis that the means of e~ is not significantly different from zero,
i

where α is the error probability,  0 is the square-root of the variance component, ̂ 0 is

the square-root of the estimated variance component, and ̂ 0i is the square-root of the
estimated variance component excluding the suspected outlier in the i-th observation, and
r denotes the degree of freedom. When a large data sample (i.e., over 100 observations) is
used, equation (4.24) and (4.25), together with their critical values, can be approximated
by equation (4.23) (Kavouras, 1982, pp.44–48).
Extra large data samples (i.e., over 1000 observations) are common in satellite
geodesy measurements nowadays. It renders the equations (4.23) through (4.26), more or
less, impractical because Q~e is computationally expensive for the matrix multiplication
when the number of observations becomes too large as can be seen from equation (4.5).
Since e~ is an underestimate of e, so, too, Q~e is an underestimate of Qe  Q y as
manifested from equations (4.4) and (4.5).
Therefore, the following modification of the Thompsonized residual of equation
(4.24), despite being an approximation, is applied in common practice (van Loon, 2008)
via
e~i n
vi    i   vi (4.27)
ˆ 0 (Q y ) ii r
because ̂ 0 is an indicator of goodness-of-fit in its squared form instead of  0 which is
usually unknown. Note that this statistic has a slightly different distribution than the
Pope’s τ distribution, since ̂ 0 and Q y are used. The sensitivity is also less than that of
equation (4.23) to (4.25), because Q~e is smaller than Q y , in the Lӧwner partial ordering;
as manifested from equation (4.5). However, the computations are extremely fast since
all quantities have been computed beforehand. This test statistic is used for the outlier
57
detection and the down-weighting incorporated into the robust estimation technique
described below.

Robust estimation
The M-estimator was introduced by Huber (1964) as a generalized form of maximum
likelihood estimator. Consider a sample of independently distributed random
variables y1 , y 2 ,..., y n from the probability density function f ( yi  Ai  ) , where  is the
vector of location parameters that are sought, with ei  yi  Ai  . The maximum
likelihood function for the estimation of  is replaced by the minimum of the log-
likelihood as
n
M  ln L ( )    (ei )  min (4.28)
i 1

where  (ei )   ln f ( yi  Ai  ) . Hence, differentiating M with respect to the location

parameters  gives subsequently
M n
 (ei ) ei n
 ( y i  Ai  ) n
   (ei )    Ai  (ei )
T
 (4.29)
 i 1 ei  i 1  i 1

and
~
T  ( ei ) ~
n n

  Ai w(e~i )~
ei  0 for e~i : yi  Aiˆ
T
Ai ~ ei  (4.29a)
i 1 ei i 1

in order to minimize M. Equation (4.26a) can be rewritten in matrix notation, in which

the normal equations are obtained as
AT P (~  
e )e~  A T  P (e~ )  y  Aˆ  0 (4.30)
such that the vector of the location or unknown parameters is found as

ˆ  AT P (~ e ) A AT P (~
1
e)y for e~ : y  Aˆ. (4.31)
~
Because the weight matrix P(e ) , in general, depends on the residuals which are
unknown, this equation may have to be solved through iterative reweighting as follows:
 
ˆ ( k 1)  AT P (e~ ( k ) ) A AT P (e~ ( k ) ) y, for k  0,1,2,...
1
(4.32)
e~ ( k 1)  y  Aˆ ( k 1) (4.33)
~
e ( k 1) P (e~ ( k ) )e~ ( k 1)
T
2 ( k 1)
ˆ0  (4.34)
nm
in the (k+1)-th iteration, after a properly chosen matrix P (~ e (k ) ) that is responsible for
down-weighting the observations according to the size of the normalized value vi . The
computation is initially started with P (~ e ( k ) )  P when k = 0. A new weight matrix is
then computed from the residuals using equation (4.35) below iteratively. The iteration is
( k 1) (k )
repeated until the difference ˆ 0  ˆ 02   , in addition to the conventional practice
2

ˆ ( k 1)  ˆ ( k )   , where δ is a small threshold value, because equation (4.31) is a global

58
indicator of goodness-of-fit particularly for a large number of observations. Note also that
the parameters with a small value in the vector ˆ can be very sensitive to changes in the
weights, potentially causing a problem for the convergence of the overall vector ˆ .
Previous comparative studies on different weight updates showed a mixed
performance (Hekimoglu and Erenoglu, 2007; Knight and Wang, 2009; Sisman, 2010),
except for the weight updates developed by Yang (1999) that exhibited an average
performance amongst all existing weight functions for the robust M-estimation. This
weight update, together with equation (4.24), is given as follows:
 P (e~ ( k 1) ) , if vi  c 0
 2
 c c v 
P (e~ ( k ) )   P (e~ ( k 1) )  0   1 i  , if c 0  vi  c1 (4.35)
 vi  c1  c 0 
 0 , if vi  c1

where c0 and c1 are chosen as 2.57 and 4.00 respectively in this study.
The above described methodology is employed for global ocean tide modeling, using
multi-mission altimeter data with distinct spatial and temporal resolution and different
tidal aliasing characteristics. We will describe the multi-altimetry data used and the
computational procedures involved in the next subsection.

4.2.5 Altimetry Data Used and Computational Procedures

The tidal analysis is conducted using multi-mission radar altimeter data from TOPEX,
Jason-1, GFO, and Envisat, which have different inclinations, and hence, different spatial
and temporal resolutions (Table 4.3). Jason-2 data is not included in the tidal analysis but
it is served for an external validation of the ocean tide models through sea surface height
(SSH) anomaly variance reduction test as will be described in Chapter 5.
The ground tracks of altimetric satellites do not repeat exactly every cycle but they
drift within 3-km along-track and 1-km cross-track with respect to nominal ground tracks
(Yi, 1995), due to non-gravitational perturbing forces such as air drag and solar radiation
pressure on the satellite. 1-Hz SSH has to be reduced to nominal ground track locations
via gradient correction because these gradients across the drifting ground tracks cause an
apparent increase in oceanic variability (Brenner et al., 1990; Dorandeu et al., 2003).
Preprocessing of those data was made through the updates and the retrieval of the
stackfile data system. This system provides gradient-corrected and edited data at their
respective nominal locations (Guman, 1997; Kruizinga 1997; Yi, 2010) of the
Geophysical Data Record (GDR). The average value (i.e. MSS) is subtracted from the
reduced 1-Hz SSH data time series at each location along the satellite altimetry ground
track. Note that the MSS of each location is updated whenever new cycles of data
become available. Mission specific correction models and better orbits are upgraded
when they become available.

59
 Table 4.3 Altimeter data used for the ocean tide modeling

Mission/Phase Time Span Cycles Source

TOPEX/POSEIDON Aug 1992 – Dec 2005 4 – 479 MGDR-B (NASA)
GFO Jan 2000 – Nov 2007 37 – 204 GDR NOAA
Jason-1 (version C) Nov 2001 – Apr 2010 1 – 303 GDR-C (CNES/NASA)
Envisat Mar 2002 – Jul 2009 10 – 80 GDR ESA/CNES
Jason-2 June 2008 – Apr 2010 1 – 66 GDR-C (CNES/NASA)

A geophysical model of oceanic response to high-frequency atmospheric variations,

called Dynamic Atmospheric Corrections (DAC) (Carrère and Lyard, 2003; Pascual et
al., 2008), is used to replace the static inverted barometer (IB) correction. The deviation
of DAC from the static IB correction is apparent along most coastal regions or in-land
seas (Figure 4.4). This indicates that a potential improvement of SSH anomaly data is
anticipated, provided that DAC is indeed a more accurate representation of atmospheric
pressure forcing.
In addition to other standard corrections, the SSH anomaly data time series is also
corrected with the ocean tidal loading effect using the NAO.99b ocean tidal loading
model, the radiational potential effect (Cartwright and Ray, 1994), and the free core
nutation (FCN) resonance effect (Matsumoto et al., 2000), as described in section 4.1.2.

Figure 4.4 The standard deviation of the difference between Static Inverted Barometer (IB) Correction and
Dynamic Atmospheric Correction (DAC)

Because the tidal aliasing effect should be mitigated by a shorter effective sampling
period, altimeter data from different satellite nearby adjacent and crossing ground tracks
60
are acquired at each grid center through a search area of 0.75o×0.75o around the grid
center for the tidal analysis (Figure 4.1). The search area at this size preserves locality of
ocean tides, since ocean tides can be very different even for a few kilometers apart,
particularly near the coasts. This additional sampling from nearby adjacent and crossing
ground tracks provides useful information for a better tidal signal recovery (Knudsen,
1994; Smith, 1999). Note that the neighboring grids may use partly the same data due to
the search area for each grid center separated by 0.25o as shown in Figure 4.5. This
partial overlapping use of data among the neighboring grid points generates the spatially
correlated tidal estimates among grid points, and hence, ensuring the coherence and
smoothness of the tidal estimates, because tides should, by nature, be smoothed. Any
abrupt spatial change should be attributed to improper handling of biases among its
neighborhood, such as satellite biases. In addition, no solution is attempted at those grid
centers where insufficient data are present (i.e., number of observations is less than three
times the number of parameters) or only Envisat data are available because of fewer data
cycles and its aliasing effect, as discussed in section 3.2. For computational efficiency
and simplicity, spatial and temporal correlations among data are not empirically
estimated, though they can be obtained as mentioned from the beginning of section 4.2.2.
Grids without tidal estimates are filled through interpolation. For higher latitude
regions (i.e.,   66 ), they are filled by the GOT4.7 model. In the process, General


Mapping Tool (GMT) land mask at 1/16o×1/16o gridded resolution was utilized to
identify the land surface so that no tidal estimates are filled onto the land surface. Hence,
an initial ocean tide model of this dissertation is determined.

Figure 4.5 Concept of the tidal solution approach

61
 As discussed in section 3.1.2, though the observed SSH anomaly contains general
ocean circulation signals, tidal signals, and errors (geophysical and media corrections)
and data noise, the observation model assumes the remaining error as random noise. This
causes the non-linearity in the solution process. To mitigate the non-linearity in the
estimation procedure, an incremental tidal analysis is conducted to estimate the remaining
tidal signals from the initial determined (background) ocean tide model using the same
combination technique, More explicitly, the incremental observations,
hssha ( ,  , t )   pred ( ,  , t ) , are utilized for solving incremental corrections of ocean tides
obs

at predefined gridded location (  c , c ), where  pred ( ,  , t ) is the tidal height predicted

from initial gridded ocean tide model by bilinear interpolating the tidal constants to the
actual observed location (  ,  ) along the satellite altimeter ground tracks. The effect of
initial and incremental tidal analysis can be seen in Figure 5.6 and Figure 5.7.
Note also that TOPEX, GFO, and Envisat ground tracks are different in locations as
well because of the orbit design. It should be pointed out that the SSH anomaly
measurements have been used twice, i.e. once in the initial ocean tide model, and once in
the incremental tidal analysis. This should be avoided in a typical Least-squares solution
(Chen, 2007, pp.33-34), so we can only apply the procedure above when a better tidal
solution can be achieved.
Standard deviation of tidal constant estimates is a good indicator to decide if the
follow-up incremental tidal analysis is of good quality. The incremental tidal constants
with standard deviation larger than ±10 cm were not included in the final solution. Any
abrupt changes in magnitude of incremental tidal constants were excluded for the final
solution.

4.3 Chapter Summary

We provide a comprehensive review on the global ocean tide modeling approaches
using satellite altimetry. We focus on formulating the problem for the observation
equation setup and the novel spatio-temporal combination approach for weighting
multisatellite altimetry data both in space and in time simultaneously based on the
specification of spatial and temporal (co-)variances in the weighted least-squares solution
process. We also discuss the possibility to include an empirical block-diagonal structure
of the cofactor matrices for both the spatial and the temporal covariances. The robust
estimation technique is also utilized for potential outlier down-weighting in the least-
squares solution process. The complete computational procedures, the usage of multi-
satellite altimetry data, and their corrections, are detailed.

62
 Chapter 5: Ocean Tide Modeling Solutions and Accuracy Analysis
In this chapter, we show the result of two global ocean tide models based on the
spatio-temporal combination using multiple mission satellite radar altimetry data.
Particular emphasis is paid on the accuracy assessment both globally and regionally,
when compared to other contemporary models.
The sea ice has been melting rapidly in the Arctic Ocean during the past two decades.
While most research effort is paid on the climate variability and the sea-level trend
projection due to sea ice melt, accurate tidal prediction in the ice-covered polar oceans
remains elusive. Potential evidences for seasonality of ocean tides based on altimetry data
is also studied in this chapter.

5.1 External Accuracy Assessment Using Tide Gauges and Tidal Solutions under
Different Weighting Schemes
To assess the accuracy of ocean tides models, tide-gauge-derived harmonic constants
(i.e. amplitude and phase or equivalently in-phase and quadrature amplitudes defined in
section 2.3.2) are necessary to serve as a reference ground truth since their longer
temporal sampling when compared to the altimetry-derived harmonic constants. Two sets
of ‘ground truth’ harmonic constants, provided by Dr. Richard D. Ray, were used for this
purpose: (1) Deep ocean tidal constants derived from 49 island and 53 bottom pressure
recorder (pelagic) stations, all together 102 sites 9 (Shum et al., 1997); and (2) coastal
ocean tidal constants derived from 739 coastal sites (Figure 5.1).

9
In ocean tide community, the harmonic constants derived from these 102 sites situated in deep ocean,
which refers to ST102p dataset, are called ‘pelagic tidal constants’.
63
 Figure 5.1 The locations of both pelagic (white stars) and coastal (red triangles) tide gauges.

Since the final gridded tidal constants vary in space only, the final tidal harmonic
constants near the tide gauges (i.e. ground truth values) can be bilinear interpolated for
the evaluation of the gridded ocean tide models. Note that not all eight major tidal
harmonic constants are present in the provided tide gauge records, especially near the
coasts.
The evaluation was made by computing the RMS difference of harmonic constants
for each constituent k generated from an ocean tide model against the reference ground
truth data, which is defined as (Andersen, 1995b):

RMS k 
1 N

2 N i 1
   2

C ksol (i)  C kref (i)  S ksol (i)  S kref (i)
2
(5.1)

where C ksol (i ) , C kref (i ) , S ksol (i ) and S kref (i) are the in-phase and quadrature amplitudes
terms, defined in equation (2.34), for the ocean tide solution and the reference ground
truth data respectively for each location i. N is the total number of locations where the in-
phase and quadrature amplitudes are computed. Root Sum of Squares (RSS), accounting
for the total deviation of the M major constituents for each model against the reference
ground truth data, is an indicator of the overall discrepancy of the model against the
reference ground truth, which is defined as:
M
RSS   RMS
j 1
2
j (5.2)

The Root Sum of Squares of the In-phase and Quadrature amplitudes (RSSIQ) for the
reference ground truth data over M major constituents is also computed. This served as a
denominator for the assessment of the overall fraction of error of the ocean tide models
against the ground truth data obtained from RSS, which is defined as:
64
 RSSIQ 
1 M N
 
2 N j 1 i 1
  2

Ckref (i)  S kref (i)
2
(5.3)

As a consequence, discrepancy D is defined as a relative error between the model and the
tide-gauge derived harmonic constants, which can be computed as:
RSS
D  100 % (5.4)
RSSIQ
Larger values of D indicate larger error in the tested ocean tide models against the ground
truth data.
This section is divided into two parts. First, the aforementioned tide gauges are used
to test the effectiveness of the proposed spatio-temporal combination approach against
the equally weighted solution, the weighted solution based on spatial (co-)variances, and
the weighted solution based on temporal (co-)variances under the same settings. Only
then, the global ocean tide models are made since it is a computationally demanding task
as mentioned in section 4.2.2. Second, we provide an accuracy assessment of two newly-
made ocean tide models (i.e. OSU12sw and OSU12vce) as compared against
contemporary ocean tide models using these reference ground truth data.
Recall that the spatio-temporal combination approach employed in this study requires
the value  to be empirically determined for the best scale of correlation length.
Therefore, ocean tidal gridded solutions near coastal tide gauges are estimated using the
proposed approach and spatial weight based on specified (co-)variances. The reason of
using coastal tide gauges only is that the proposed approach has a negligible effect in
Deep Ocean as shown in Figure 4.3. It was found that the best  values for OSU12sw,
OSU12vce, and the weighted solution based on spatial (co-)variances are 90, 40, and 30,
respectively, based on the result from initial and incremental tidal analysis (Figure 5.2).

65
 Figure 5.2 Ground truth comparison of ocean tide solutions near 507 coastal tide gauges based on initial
(upper) and incremental (lower) tidal analysis using spatio-temporal combination approach(i.e. OSU12sw
and OSU12vce) and the weighted solution based on spatial (co-)variances (i.e. spatial weight).

66
 Figure 5.3 Ground truth comparison of ocean tide solutions near 102 pelagic and 507 coastal tide gauges
for the spatio-temporal combination approach(i.e. OSU12sw and OSU12vce), the equally weighted
solution (i.e. equal weight), the weighted solution based on spatial (co-)variances (i.e. spatial weight), and
the weighted solution based on temporal (co-)variances (i.e. temporal weight) based on initial and
incremental tidal analysis.

These respective  values are subsequently employed for generating the global
gridded ocean tide models and testing against the equally weighted solution, the weighted
solution based on spatial (co-)variances, and the weighted solution based on temporal
(co-)variances. It was found that spatio-temporal combination approach has a good
performance in both Deep and Shallow Ocean among different weighting schemes
(Figure 5.3).
Note that OSU12sw model is based on conjectured weight component estimates for
the covariance matrix, the test of the OSU12sw solution robustness has also been
conducted to examine the changes using the coastal tide gauges. If the change is
negligible, the solution is robust. This test is conducted by distorting the estimated weight
components, ̂ S , ̂T , and ̂ ST , by ±5%, ±7%, and ±10%, individually. It was found that
the change of the solution is negligible, no matter for initial and incremental tidal analysis
(Figure 5.4).

67
Figure 5.4 Robustness test for the OSU12sw solutions using 507 coastal tide gauges under the distortion of
~ S2 (Var_S), ~T2 (Var_T), and ~ST2 (Var_ST), by ±5%, ±7%, and ±10% for initial and incremental tidal
analysis.

Eleven contemporary ocean tide models were included in the comparison: DTU10
(Cheng and Andersen, 2011), EOT08a/10a/11a (Savcenko and Bosch, 2008), FES2004
(Lyard et al., 2006), GOT00.2/4.7 (Ray, 1999), TPXO6.2/7.1/7.2 (Egbert and Erofeeva,
2002), OSU12sw and OSU12vce (this study). Pelagic and coastal tide gauges were used
where bilinear interpolation of the model harmonic constants is possible for all current
ocean tide models.
It should be noted again that some of the ocean tide models included tide gauge (and
altimetry) data either as constraints (or assimilated into hydrodynamic models to compute
ocean tides) or directly used the data to estimate tides. Thus the evaluation using ground-
truth tide gauges may not be completely independent for some of the models, e.g.,
FES2004 or the TPXO models. In addition, the empirical models (i.e. EOT08a/10a/11a,
GOT00.2/4.7, and DTU10), which used these models for incremental tidal analysis,
would have much higher spatial resolution and better coverage in coastal regions, as
FES2004 model affords higher spatial resolutions than purely empirical ocean tide
models.
The two purely empirical models in this study, OSU12sw and OSU12vce,
demonstrate a comparable performance in terms of accuracy assessment when compared
to other ocean tide models using pelagic tide gauges. M2 tide is shown to have a
68
substantial improvement (Table A.1 in Appendix A). However, the accuracy is still at 2–
3 cm level, representing ~7.5% relative error with respect to the total tidal signal of the
pelagic tide gauges. This is primarily due to tidal aliasing effect into the ocean tidal
solution arising from under sampling. When it is compared to other models using coastal
tide gauges, the models in this study have a substantial improvement over other models
(Figure 5.5), which indicates an improvement of ~2–3% over GOT4.7 model (i.e. second
best model). Note that OSU12vce has ~1500 gridded estimates (i.e. <0.3%) less than
OSU12sw which is due to numerical instability arising from the variance component
estimation (VCE) algorithm. OSU12sw is also slightly better than OSU12vce.

Figure 5.5 Ground truth comparison of ocean tide models at both 102 pelagic and 549 coastal tide gauges.

5.2 Internal Accuracy Assessment Using Variance Reduction Test

Tide is the major signal in ocean. The best model should minimize most variation of
SSH anomaly residuals (King et al., 1995). Therefore, standard deviation of SSH
anomaly along satellite tracks and that of SSH anomaly residuals after the removal of
tidal height predictions of ocean tide models are computed over the entire sea region to
investigate how much the oceanographic variability be minimized in this assessment.
The SSH anomaly residual of location (,  ) at time epoch time t, hssha res
( ,  , t ) , is
defined as:
res
hssha ( ,  , t ) = hssha
obs
( ,  , t ) – (diurnal + semidiurnal tides) – LP (5.5)

69
 obs
where hssha is the observed instantaneous SSH anomaly as defined in equation (3.3), the
diurnal and semidiurnal tides are predicted by the ocean tide models, and the equilibrium
long-period (LP) tides are calculated based on Cartwright and Edden (1973) (which was
adopted in the GOT00.2 and GOT4.7 models) for consistency.
Note that the gridded global ocean tidal constants of diurnal and semi-diurnal
constituents are bilinear interpolated to the positions of SSH anomaly along satellite
obs
tracks before the prediction and subtraction from hssha . Hence, the variance of SSH
obs res
anomaly hssha and its residual hssha after ocean tide correction for each location s along
satellite altimeter tracks are computed as
ns

 (h
j 1
obs
ssha ( s ,  s , t j )  hssha
obs
( s ,  s )) 2
obs
var(hssha ( s ,  s , t j ))  (5.6)
ns  1
and
ns

 (h
j 1
res
ssha ( s ,  s , t j )  hssha
res
( s ,  s )) 2
res
var(hssha ( s ,  s , t j ))  (5.7)
ns  1
where hsshaobs
( s ,  s ) and hssha
res
( s ,  s ) are the mean SSH anomaly and its residual
respectively, and ns is the number of observations at location s along satellite altimeter
obs res
ground tracks. Consequently, the standard deviation of hssha and hssha within the entire sea
region of each altimeter p can be computed, respectively, as
mp ns

 (h obs
ssha ( s ,  s , t j )  hssha
obs
( s ,  s )) 2
s 
obs
ssha p 
s 1 j 1
mp
(5.8)
 (n
s 1
s  1)

and
mp ns

 (h res
ssha ( s , s , t j )  hssha
res
( s , s )) 2
s 
res
ssha p 
s 1 j 1
mp
(5.9)
 (n
s 1
s  1)

where mp is the total number of along satellite ground track locations for each altimeter p
within the entire sea region. The SSH variance explained (VE) by ocean tides for each
altimeter data is computed as:
s ssha
obs
p  s ssha
res
p
VE p 
p  100%
(5.10)
s ssha
obs

Hence, the overall SSH variance explained (VE) by ocean tides for all altimeter data
are computed by:

70
  s   s 
q q
obs 2 res 2
ssha p  ssha p
p 1 p 1
VE   100% (5.11)
 s 
q
obs 2
ssha p
p 1

where q is the total number of satellite altimeters.

Figure 5.6 Overall sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide
models. Note that deep and shallow oceans are bathymetry depth greater than 1000 m and less than 1000
m, respectively. Jason-2 altimeter data serves as an independent data for validation.

The overall statistics of the deep, shallow, and overall ocean in the variance reduction
test was summarized in Appendix A (Table A.3, Table A.4, and Table A.5). Figure 5.6
displays the SSH variance explained (VE) by the ocean tides for the deep, shallow, and
overall ocean. It indicates OSU12 models in this study minimize the most variation of
SSH anomaly and hence, attaining the maximum variance explained by ocean tide as
manifested in its percentage. Note also that Jason-2 altimeter data was excluded in the
generation of the OSU12 model solutions and therefore, serving an independent data used
for the evaluation when it is compared to DTU10 and EOT08a/10a/11a models in which
Jason-2 data has been included in the models already. Nevertheless, it is still comparable
to all other models.

71
 It is surprising that ~78% of SSH anomaly variation is explained by the ocean tides in
Shallow Ocean when compared to that of Deep Ocean (i.e. ~70%), since it is expected
the barotropic response of ocean circulation and hydrodynamic effects near shallow sea
would dominate the SSH variability. As a result, it is worthwhile investigating the
percentage of SSH variability explained in Shallow Ocean, particularly over high
hydrodynamic regions. This serves as a guideline for those who are ultimately interested
in modeling regional ocean dynamics. The result will be shown in the next sub-section.
To demonstrate the uniform improvement of the OSU12 model solutions for the
globe, we produce indexes sorted by the standard deviation of SSH anomaly residuals per
grid for visual comparison purpose besides the overall statistics. The model yielding the
minimum SSH anomaly residual variation is assigned to 1 whereas the one with the
maximum SSH anomaly residual variation is assigned to 11, because eleven models were
included in the comparison. The results are shown in Figure 5.7 (i.e. OSU12sw), Figure
5.8 (i.e. OSU12vce), Figure 5.9 (i.e. DTU10), Figure 5.10 (i.e. EOT11a), Figure 5.11 (i.e.
GOT4.7), Figure 5.12 (i.e. FES2004), and Figure 5.13 (i.e. TPXO7.2). Note that a further
reduction of SSH anomaly residual variation for OSU12sw and OSU12vce models is
achieved via incremental tidal analysis (Figure 5.7 and Figure 5.8).

72
Figure 5.7 Index of the minimum SSH anomaly residual variation (OSU12sw) before (upper) and after
(lower) incremental analysis

73
Figure 5.8 Index of the minimum SSH anomaly residual variation (OSU12vce) before (upper) and after
(lower) incremental analysis

74
Figure 5.9 Index of the minimum SSH anomaly residual variation (DTU10)

Figure 5.10 Index of the minimum SSH anomaly residual variation (EOT11a)

75
Figure 5.11 Index of the minimum SSH anomaly residual variation (GOT4.7)

Figure 5.12 Index of the minimum SSH anomaly residual variation (FES2004)

76
 Figure 5.13 Index of the minimum SSH anomaly residual variation (TPXO7.2)

5.3 Global Result and Coastal Assessment

No apparent difference is observed for OSU12sw and OSU12vce models, while the
difference of M2 tide between OSU12sw/vce and EOT11a ocean tide model as shown in
Figure 5.14. The Greenwich phase lag contours at coastal regions for EOT11a model are
denser than that of OSU12sw/vce, due to finer spatial resolution of EOT11a model (with
FES2004 as background model). Noticeable difference in the pattern of amplitude is
found, e.g. Indian Ocean.

77
 Figure 5.14 Global M2 co-tidal chart of OSU12vce (upper) and EOT11a (lower) ocean tide models.

Besides the global result, coastal or shallow sea regions with high dynamic oceanic
variability are of great interest to physical oceanographers. A regional assessment was
also conducted for the ocean tide model evaluation study.
Six coastal or shallow sea regions with high dynamic oceanic variability were chosen
in the tide model evaluation study. They are the Gulf of Mexico and Northwest Atlantic
(GMNA), Patagonia Shelf, Southeast Australia, Indonesia, Northeast Pacific, and Sea of
78
Japan coastal regions. These regions are shown in Figures 5.15 through 5.20, with the
standard deviation of multisatellite altimetry SSH anomaly residuals using OSU12vce
model as the tide correction, together with locations of tide gauge sites with pelagic tidal
constants used for tide model comparisons.
The percentage for the discrepancy D of the ocean tide models against the selected
ground truth data, based on RSS and RSSIQ, were calculated to assess the fraction of
error. Comparison of this percentage among the ocean tide models reveal a large
disagreement between the ocean tide models and the tidal records at coastal sites,
particularly in the GMNA, Northeast Pacific and Southeast Australia regions, where the
disagreement exceeds 40% (Figure 5.21). The disagreement of other regions is, on the
other hand, ~25% in general. This implies the tidal variability in shallow water is not well
represented.
The ocean tide models, when compared to the coastal tide gauge constants, display
heterogeneous performance over Patagonia Shelf and Southeast Australia regions,
showing different approaches in handling the regional hydrodynamics near the coast, in
particular FES2004 and TPXO6.2/7.1/7.2 models since these models were generated
from hydrodynamic assimilation approach.
DTU10, EOT08a/10a/11a, and GOT4.7 models show an improved result over
FES2004 and GOT00.2 model, since DTU10 and EOT08a/10a/11a models are indeed
based on FES2004 as a-priori model and GOT4.7 model is a successor of GOT00.2.
TPXO6.2/7.1/7.2 models exhibit heterogeneous performance when compared with tidal
constants in coastal sites, depending on the investigated regions. Overall, the
OSU12sw/vce models demonstrate an ‘above average’ performance when compared to
other ocean tide models as against the coastal tide gauge constants, in particular
outperform in Southeast Australia and Sea of Japan regions.
Comparison of tidal constants of these models with in situ data at respective pelagic
sites (depicted with star symbol in the Figures 5.15 – 5.20) was conducted in a similar
fashion (Figure 5.21). The percentage of the disagreement among ocean tide models
reveals a relatively better agreement between the ocean tide models and the tidal constant
at pelagic sites than that at coastal sites, where the disagreement is less than 10% in
general except Indonesian sea and Sea of Japan. The ocean tide models, compared
against the pelagic tidal constants, demonstrate homogeneous performance for all the
regions selected in this investigation, except for Sea of Japan regions. It is also
indistinguishable which version of TPXO models could provide better result when
compared with tidal records at pelagic sites, depending on investigated regions.
Note that pelagic sites are scarce in all the study regions. For example there are only 2
pelagic sites around Sea of Japan. As a result the analysis associated with pelagic data
test in these sites could be statistically insignificant.

79
Figure 5.15 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Gulf of Mexico and Northwest Atlantic region.

Figure 5.16 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Patagonia Shelf region.

80
Figure 5.17 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Southeast Australia region.

Figure 5.18 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Indonesian Sea region.

81
Figure 5.19 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Northwest Pacific region.

Figure 5.20 Standard deviations of SSH anomaly residuals at along satellite track locations and location of
tide gauge sites used as the ground truth (pelagic sites shown as stars and coastal sites as triangles) in the
Sea of Japan region.

82
Figure 5.21 Ground truth comparison of ocean tide models at pelagic (upper) and coastal (lower) tide
gauges with number of tide gauges indicated in x-axis.

83
Figure 5.22 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Northeast Pacific region. Jason-2 altimeter data serves as an independent data for validation.

Figure 5.23 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Gulf of Mexico and Northwest Atlantic region.

84
Figure 5.24 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Patagonia Shelf region.

Figure 5.25 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Sea of Japan region.

85
Figure 5.26 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Indonesian sea region.

Figure 5.27 Sea surface height (SSH) variance explained by ocean tides (in percent) for ocean tide models
in Southeast Australia region.

86
 It was shown in Figure 5.6 that the SSH variability explained by ocean tides for
shallow water regions is ~80%. However, Sea of Japan where ocean tide models present
the least SSH variance explained by ocean tides (i.e. ~51%). As a matter of fact, all
regions show heterogeneous variability explained by ocean tide corrections for the deep
ocean SSH variance reduction study (Figure 5.22 through Figure 5.27). This is mostly
because the areas with high standard deviations of SSH anomaly residuals are located in
the western boundary current (e.g., Gulf Stream) (Figure 5.15). For example, in GMNA
region, the main cause of this result is the transport of warm water from Caribbean Sea
through Yucatan Channel generates the loop current in the eastern Gulf (east of the ca.
272.5oE longitude) (Sheinbaum et al., 2002). This loop current eventually spreads and
forms anticyclonic (warm-core) eddies at the central of the Gulf (ca. 266 – 272.5oE
longitude) and their associated cyclonic (cold-core) eddies (ca. west of 266oE longitude);
these are the primary circulatory features of the region (Davis et al., 2001) as could be
seen from Figure 5.15 of this study, and from Plate 3 of Leben et al. (1990) that used
Geosat altimetry data without ocean tide correction. As a consequence, the interaction of
ocean tides with the aforementioned non-tidal circulation features causes tidal mixing due
to the shape of this ocean basin. Given the coarse temporal resolution of satellite
altimetry, the altimetry observations are unable to capture or separate the tides accurately,
not to mention the tidal mixing phenomena.
Overall, the OSU12sw/vce model shows a minimum residual oceanographic
variability when compared to other models, both in Shallow and Deep Ocean.

5.4 Seasonal Variations of Ocean Tides in the Sub-Arctic Ocean

Sea ice in the Arctic Ocean is sensitive to small changes in vertical oceanic heat flux,
and is undergoing rapid thinning during the past decades (Holloway and Sou, 2002;
Kwok and Untersteiner, 2011). Its spatial coverage and change or diminishment in Arctic
Ocean is important information for assessing large-scale climate changes (Renganathan,
2010). The measurement of the sea-ice freeboard height change and its mapping via
satellite laser or radar altimetry is the quantity for describing its thickness and coverage.
However, accurate knowledge of dynamic ocean topography, the geoid, and ocean tides
are required (Kwok and Untersteiner, 2011; Renganathan, 2010).
Indeed, there are evidences that seasonally sea-ice-covered ocean causes seasonality
of tidal signals (Kagan et al., 2008; Renganathan, 2010), which would affect the accuracy
of sea-ice freeboard height change retrieval using satellite altimetry. The sea-ice cover
also has a substantial effect on ocean tides in the Arctic Ocean. Recent simulation study
on ice-induced changes of the tidal dynamics over the Siberian continental shelf shows
that the sea-ice covers caused a decrease in tidal amplitudes and an increase in tidal
phases, and hence changing the tidal dynamic characteristics and its energy dissipation
(Kagan et al., 2008). A more recent simulation study on seasonal variations of M2 tide in
the Arctic Ocean concluded that the change in tidal amplitude and in phase do not exceed
a few centimeters and a ten of degree (i.e., 10o) for summer and winter seasons,
respectively (Kagan et al., 2011).

87
 In addition to sea-ice, significant seasonal modulation of ocean tides was found over
high latitude region (e.g. European shelf), as revealed by using TOPEX altimeter and tide
gauge data to estimate the seasonal modulation constituents of ocean tides
(Leeuwenburgh et al., 1999; Huess and Andersen, 2001). However, the aforementioned
studies focus on either the use of tidal hydrodynamic models to either simulate or
examine seasonal variations of M2 tide or seasonal modulation of M2 tidal constituent in a
European Shelf where strong tidal mixing and surge are present. No experiment study has
been conducted to reveal the potential seasonality of ocean tides using satellite altimetry
data.
Here we use multi-mission satellite altimetry at latitude between 60o and 72o in the
sub-Arctic Ocean, with the data span tabulated in Table 4.3, to test the hypothesis
whether plausible evidences of the seasonality of tidal signals exist in the seasonally ice-
covered sub-Arctic Ocean through empirical ocean tide modeling.
Indeed, the feasibility of an empirical ocean tide modeling solution using GFO and
Envisat around Antarctica Ocean has been demonstrated (Yi et al., 2006). It was shown
to have a comparable accuracy with the TPXO6.2 and FES2004 models, which are
hydrodynamic models with data assimilation at latitude higher than ±66o. Similarly, it can
be applied in Arctic Ocean. The altimeter data used and preprocessing for this study can
be found in section 4.2.5.
To quantify and detect the seasonal variation of ocean tides, tidal analysis
methodologies described in Chapter 4 and variance reduction test (i.e., internal accuracy
assessment) are conducted, with the implicit assumption that the variance reduction
ability by ocean tide models is different during summer and winter season. Notice that we
define the winter season between September and March, whereas summer season is
defined between March and September per year for the separation of data time span. The
data amount per grid point is sufficient for robustness of the seasonal ocean tidal solution
owing to the altimeter tracks converging in the polar ocean.
Seven global and regional ocean tide models were used for the aforementioned
purpose: DTU10 (Cheng and Andersen 2011), EOT11a (Savcenko and Bosch 2008),
FES2004 (Lyard et al. 2006), GOT4.7 (Ray 1999), NAO.99b (Matsumoto et al. 2000),
TPXO7.2 (Egbert and Erofeeva 2002), and the AOTIM-5 regional tide model (Padman
and Erofeeva 2004). All the tide models did not use Jason-2 altimetry data in the solution,
except for the DTU10 and the EOT11a models. In other words, Jason-2 data are not
independent in the evaluations of these two tide models.
All (6) global ocean tide models exhibit similar performance, except for AOTIM-5
regional tide model, which has a substantially worse performance (Shum et al. 2006)
(also Table 5.1). Multi-mission altimetry SSH anomaly residuals after tidal corrections
using the 7 global and regional tide models is at 9–12 cm RMS (the regional model
AOTIM-5 is excluded in the residuals computation) in the Subarctic Ocean. The overall
SSH anomaly variation (before and after ocean tide correction) shows substantial
differences between summer and winter seasons (Tables 5.2 and Table 5.3), with larger
SSH anomaly variation for all altimeter data during winter season when compared to that
of summer. Bold values in Table 5.1 through Table 5.3 are associated with the tide
models with the best performance. The variation of SSH anomaly residuals during
88
summer season is lower than that of during winter season by 15–30%, no matter it is
before or after ocean tide correction. This is particularly evident in TOPEX and Jason-1
interleave, and Jason-2 altimeter data, since the data time span covered less than three
years only.
A variance reduction test was also conducted for the seasonal tide models (Tables 5.2
and 5.3). Except TOPEX and Jason-2 altimeter data primarily due to ±66o coverage, our
investigation focuses on a higher latitude region. Substantial improvement in the variance
reduction of all altimeter data during summer season was made when compared to the
empirical model without seasonal tides determination in this study (Table 5.2 and Table
5.3). Only GFO and Envisat demonstrate an improvement during winter probably due to
the fact that major tidal variation occurred around Chukchi Sea where latitude is higher
than 66o, where it is beyond the range of TOPEX/Jason coverage. Hence, the GFO and
Envisat altimeter observations are more sensitive to the changes.
The seasonality of ocean tides is also apparent from the Jason-2 along track SSH
anomaly residual variation, particularly in Siberian continental shelf (e.g. Sea of Okhotsk
and Bering Sea) and northern Atlantic Ocean (Figure 5.28).

89
 Table 5.1 Standard deviations of SSH anomaly residuals of the empirical seasonal ocean tide model along satellite tracks in subarctic ocean (in cm) for
entire data without separation into summer and winter seasons. Stdev (before) and Stdev (after) are the standard deviation of the SSH anomaly before and
after ocean tide correction, respectively; RSS Stdev (before) and RSS Stdev (after) are the Root Sum of Squares of the standard deviation of the SSH
anomaly before and after ocean tide correction for the entire altimeter data, respectively; VE is the variance explained by ocean tide correction for the entire
region.

Stdev (after) RSS

Model TOPEX/ TOPEX Jason-1 Stdev VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
55.67 53.56 47.53 48.54 52.82 54.92 55.91
(before)
DTU10* 10.23 9.14 10.12 9.84 11.93 9.28 11.21 27.23 80.51
EOT11a* 10.33 9.15 10.30 9.86 11.91 9.32 11.24 27.36 80.41
FES2004 10.42 9.26 10.50 10.04 12.02 9.48 11.51 27.78 80.11
GOT4.7 10.27 9.14 10.10 9.85 11.90 9.27 11.24 27.24 80.50
90

NAO.99b 10.58 9.33 10.96 10.33 11.89 9.43 11.43 28.05 79.92
TPXO.7.2 10.30 9.22 10.53 10.13 11.89 9.32 11.27 27.57 80.27
AOTIM5 26.47 26.33 24.96 24.45 27.65 26.31 26.27 - -
OSU12 10.43 8.95 9.94 9.69 12.01 9.17 11.28 27.15 80.56
RSS Stdev (before) = 139.69 cm
*Model used Jason-2 data for ocean tide solutions. Bold values indicate the lowest residual and better overall variance reduction for the respective test.
 Table 5.2 Standard deviations of SSH anomaly residuals of ocean tide models along satellite tracks in subarctic ocean (in cm) during summer season.

Stdev (after)
RSS Stdev
Model TOPEX/ TOPEX Jason-1 VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
55.37 53.06 47.14 48.19 51.66 54.88 55.13
(before)
DTU10* 9.32 8.54 9.08 9.14 9.96 8.58 9.78 24.38 82.38
EOT11a* 9.45 8.51 9.34 9.16 9.99 8.61 9.86 24.57 82.24
FES2004 9.56 8.63 9.57 9.39 10.15 8.81 10.25 25.13 81.84
GOT4.7 9.37 8.49 9.00 9.14 9.95 8.53 9.82 24.34 82.41
NAO.99b 9.68 8.71 9.94 9.67 9.78 8.69 9.89 25.12 81.85
TPXO.7.2 9.42 8.59 9.51 9.46 9.94 8.64 9.84 24.75 82.11
OSU12 9.48 8.27 9.01 9.11 10.02 8.44 9.91 24.33 82.42
This study
(Summer 9.49 8.14 8.96 9.04 10.09 8.42 10.07 24.34 82.41
Model)
91

RSS Stdev (before) = 138.37 cm

*Model used Jason-2 data for ocean tide solutions. Bold values indicate the lowest residual and better overall variance reduction for the respective test.
 Table 5.3 Standard deviations of SSH anomaly residuals of ocean tide models along satellite tracks in subarctic ocean (in cm) during winter season.

Stdev (after) RSS

Model TOPEX/ TOPEX Jason-1 Stdev VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
56.26 54.57 48.20 47.56 54.73 55.15 56.86
(before)
DTU10* 10.94 9.59 10.87 10.10 13.18 9.85 12.17 29.17 79.37
EOT11a* 11.01 9.65 10.97 10.11 13.13 9.90 12.17 29.25 79.32
FES2004 11.08 9.75 11.15 10.24 13.21 10.03 12.34 29.57 79.09
GOT4.7 10.96 9.64 10.89 10.07 13.12 9.89 12.18 29.18 79.36
NAO.99b 11.28 9.81 11.68 10.44 13.24 10.05 12.42 29.99 78.79
TPXO.7.2 10.99 9.69 11.27 10.31 13.12 9.90 12.23 29.45 79.17
OSU12 11.03 9.40 10.83 10.02 13.24 9.70 12.27 29.12 79.41
This study
(Winter 11.19 9.46 10.70 9.90 13.27 9.78 12.35 29.17 79.37
92

Model)
RSS Stdev (before) = 141.41 cm
*Model used Jason-2 data for ocean tide solutions. Bold values indicate the lowest residual and better overall variance reduction for the respective test.
 Figure 5.28 Standard deviation of residual SSH anomaly during summer season (left) and during winter
season (right) for Jason-2 altimeter data (Cycle 1–66) over Subarctic Ocean with ocean tides corrected by
GOT4.7 tide model, indicating the potential seasonality of ocean tides. Note that winter season is defined
between September and March, whereas summer season is defined between March and September.

We also examine the difference in amplitude and in phase for seasonal ocean tide
solutions estimated from altimetry data during summer season and winter season. In
general, all ocean tidal constituents respond differently in the presence of potential sea-
ice coverage. A few centimeter changes in amplitude and a ten of degree changes in
phase are observed for most regions that agree well with the findings from Kagan et al.
(2011), except around Chukchi Sea near eastern Siberia region where noticeable spatial
changes in tidal amplitude and phase are observed, particularly for S2 and K1 tides at
latitude higher than 66o (Figures 5.29 and 5.30).

93
Figure 5.29 Amplitudes of M2 (top) and S2 (middle) and K1 (bottom) tides during summer season (left) and
their difference with that during winter season (right) over the Subarctic Ocean. Note that the positive
difference indicates the amplitude in summer season is larger than that of winter.

94
Figure 5.30 Phase of M2 (top) and S2 (middle) and K2 (bottom) tides during summer season (left) and their
difference with that during winter season (right) over Subarctic Ocean.

95
 In principle, the presence of sea ice should provide a damping effect on ocean tidal
amplitudes and thus an increase in tidal phases (Kagan et al. 2008). However, the
seasonal ocean tides determined from altimeter observations reveals that the S2 and K1
tidal amplitude can be increasing during winter season when compared to that in summer.
Note that the positive difference in amplitude represents the tidal amplitude in summer
season is larger than that of winter for a correct interpretation (Figure 5.28). The changes
in tidal amplitude can be larger than 10 cm, whereas the tidal phase varies significantly
around the region. We speculate the variations in the tidal solutions or the larger observed
winter amplitudes of tides could be due to substantial amount of seasonal sea ice
freeboard movement around the region, presumably due in part to energetic ocean
dynamics, with the evidence that spatial extent of the sea ice has been observed to be
changing rapidly every month (see Figure 7.8 to 7.13, Renganathan (2010) and from
National Snow and Ice Data Center (NSIDC)).
The above result provides an evidence for seasonal tidal pattern as manifested in the
variance reduction difference of altimeter SSH residual anomaly data during summer and
winter season. The subsequent empirical seasonal (winter and summer) ocean tide model
separately derived from satellite altimetry data, yields both qualitative and quantitative
explanation on the seasonal tidal pattern, in particular around Chukchi Sea near eastern
Siberia region where the presence of large amount of sea ice freeboards does contribute
to the significant changes in M2, S2, and K1 amplitude and the large variation in phase
around the region. These results are newly found and reported in this sub-section.

5.5 Chapter Summary

The OSU12sw and OSU12vce ocean tide models developed in this study shows a
comparable or a better performance among other contemporary ocean tide models, both
in terms of external and internal accuracy assessment, particularly near the coast.
We first quantify the difference in variance reduction for multi-mission altimeter SSH
anomaly residual during summer and winter seasons, with the winter residual 15–30%
larger than that of summer. Experimental seasonal ocean tide solutions derived from
satellite altimetry revealed that the recovered winter and summer tidal constituents differ
by a few cm in amplitude and tens of degrees in phase. Relatively large seasonal tidal
patterns, in particular M2, S2, and K1 tides, were found in the Chukchi Sea near eastern
Siberia region coincident with the seasonal presence of sea-ice. This result also justifies
the necessity for an improved mean sea surface model to produce more accurate satellite
altimetry sea surface height (SSH) anomaly data, which in turn will allow advancement
of ocean tide models accounting for seasonality in the Arctic Ocean.

96
 Chapter 6: Conclusion and Future Work
Ocean tides, the main variable component in ocean surface topography, are of the
scientific and practical importance over hundreds of years, including their constraints on
Earth rotation and Earth’s ocean bottom tidal dissipation, and the meridional overturning
circulation; their corrections to various spaceborne geodetic measurements allowing more
accurate quantification of climate-change signals, and their applications to ship
navigation safety and recreational beach usage.
Current advance in satellite altimetry brings about accurate surface height
observations in a synoptic mean over the past two decades. This globally sampled data
records enables significant improvement in global ocean tide modeling. Though ocean
tides have been determined accurately to within ±(2–3) cm RMS accuracy in the deep
ocean, no contemporary tide models are equally well represented to all coastal regions or
shallow seas where ocean tidal dynamics are highly complicated due to the presence of
non-linearity in dynamic tidal motion caused by local hydrodynamic processes that is
constrained from the shape and bathymetry of the coast. Despite the gain in physical
understanding of ocean tidal dynamics near the coastal regions using hydrodynamic and
assimilation modeling, these two solution approaches depend heavily on the availability
and the quality of the bathymetry and tide gauge data for different regions. Empirical
fine-tunings of the solution is required.
In this dissertation, we first developed a novel spatio-temporal combination approach
that simultaneously considers both spatial and temporal (co-)variances for weighting
multisatellite altimetry data in a balanced manner for the gridded ocean tide solutions.
This approach were examined and proven to be effective when compared to the equally
weighted solution, the weighted solution based on spatial (co-)variances, and the
weighted solution based on temporal (co-)variances under the same settings, using
independent tidal constants from pelagic and coastal tide gauges. We speculate this
approach can be further applied to other spatio-temporally obtained datasets for better
parameter estimation, provided good spatial and temporal (co-)variance model
specifications.
Two versions of purely empirical ocean tide models, called OSU12sw and
OSU12vce, are generated by a set of empirical procedures and by variance component
estimation technique, respectively. We demonstrate comparable ocean tide model
solutions when compared to contemporary models in the Deep Ocean, whereas the model
in this study indicates a further improvement along the Shallow Ocean (i.e., ~2–3%
improvement over the GOT4.7 model) in the external assessment based on independent
tide constants from coastal gauges. The model also has a good performance in
minimizing the overall sea surface height (SSH) anomaly residual.
Because of the seasonal presence of sea ice cover, an accurate tidal prediction in the
polar oceans remains elusive. A hypothesis testing on the potential seasonality of ocean
97
tides in Subarctic Ocean was conducted. We first found a substantial difference in
variance reduction for multi-mission altimeter SSH anomaly residuals during summer
and winter seasons, with the sea surface height anomaly residuals in winter (15–30%)
larger than that of summer. This result provides an evidence for seasonal tidal pattern.
Experimental seasonal ocean tide solution was derived from satellite altimetry. It reveals
that the recovered winter and summer tidal constituents differ by a few cm in amplitude
and tens of degrees in phase, in general. Particularly large seasonal tidal patterns were
found in the Chukchi Sea near eastern Siberia region where they show coincidence with
the seasonal presence of sea ice, in particular M2, S2, and K1 tides. This result also
justifies the necessity for an improved mean sea surface (MSS) model to produce more
accurate satellite altimetry SSH anomaly data, which in turn will allow the advancement
of ocean tide models by accounting for the seasonality in both the Arctic and Antarctic
Ocean.
Further improvement of ocean tide modeling also lies in a better sea-state bias,
waveform retracking along coastal regions, and over icy and sea-ice surfaces. More
current and future altimetry missions, such as CryoSat-2, HY-2A, SARAL/Altika,
Sentinel-3, Jason-3 and ICESat-2, will further enhance both the spatial and the temporal
resolution of multi-altimeter data, particularly in higher latitude regions beyond TOPEX-
class satellites (i.e.,   66 ), and thus, the spatial resolution of ocean tide model.

98
 Bibliography

Accad, Y. and C. L. Pekeris, Solution of the tidal equations for the M2 and S2 tides in the
world oceans from a knowledge of the tidal potential alone, Phil. Trans. R. Soc., A290,
235–266, 1978.
Alcock, G. A. and D. E. Cartwright, Some experiments with orthotides, Geophys. J. Roy.
Astron. Soc., 54, 681–696, 1978.
Agnew, D. C., NLOADF: A program for computing ocean-tide loading, J. Geophys. Res.,
102(B3), 5,109–5,110, 1997.
Amante, C. and B. W. Eakins, ETOPO1 1 Arc-minute global relief model: procedures,
data sources and analysis, NOAA Technical Memorandum NGDC-24, National
Geophysical Data Centre, NESDIS, NOAA, Department of Commerce, Boulder, 2009.
Amiri-Simkooei, A. R., Noise in multivariate GPS position time-series, J. Geod., 83,
175–187, 2009.
Andersen, O. B., Ocean tides in the northern North Atlantic and adjacent seas from ERS
1 altimetry, J. Geophys. Res., 99(C11), 25,557–25,573, 1994.
Andersen, O. B., Global ocean tides from ERS-1 and TOPEX/POSEIDON altimetry, J.
Geophys. Res., 100(C12), 25,249–25,259, 1995a.
Andersen, O. B., Intercomparison of recent ocean tide models, J. Geophys. Res.,
100(C12), 25,261–25,282, 1995b.
Andersen, O. B. Shallow water tides in the northwest European shelf region from
TOPEX/POSEIDON altimetry, J. Geophys. Res., 104(C4), 7729–7741, 1999.
Andrews D., A robust method for multiple linear regression, Technometrics 16, 523–531,
1974.
Arbic, B. K., S. T. Garner, R. W. Hallberg, H. L. Simmons, The accuracy of surface
elevations in forward global barotropic and baroclinic tide models, Deep Sea Res. II, 51,
3069–3101, 2004.
Awange, J. L. and F. W. O. Aduol, An evaluation of some robust estimation techniques
in the estimation of geodetic parameters, Surv. Rev. 35 (273), 146–162, 1999.
Baarda W., A testing procedures for use in geodetic networks, Publ 2(5), Netherlands
Geodetic Commission, Delft, Netherlands, 1968.
Barbosa, S. M., M. E. Silva, and M. J. Fernandes, Time series analysis of sea-level
records: Characterising long-term variability, Non-linear time series analysis in the
geosciences, Lecture Notes in Earth Sciences, 112, 157–173, 2008.
99
Barnett, V., and T. Lewis, Outliers in Statistical Data, 3rd edn. John Wiley, Chichester,
1994.
Baselga, S., Critical limitation in use of τ test for gross error detection, J Surv Eng ASCE
133 (2), 52–55, 2007.
Bertiger, W., S. D. Desai, A. Dorsey, B. J. Haines, N Harvey, D. Kuang, A. Sibthorpe,
and J. P. Weiss, Sub-centimeter precision orbit determination with GPS for ocean
altimetry, Mar. Geod., 33(S1): 363–378, 2010.
Bosch, W., Geodetic application of satellite altimetry, International Association of
Geodesy Symposia, 126, 3–21, In: Satellite Altimetry for Geodesy, Geophysics and
Oceanography, C. Hwang, C. Shum, and J. Li, Springer, 2004.
Bosch, W., R. Savcenko, F. Flechtner, C. Dahle, T. Mayer-Gurr, D. Stammer, E.
Taguchi, and P. Canceil, Residual ocean tide signals from satellite altimetry, GRACE,
gravity fields, and hydrodynamic modeling, Geophys. J. Int. 99 (178), 1185–1192, 2009.
Box, G. E. P., and G. M. Jenkins, Time Series Analysis: Forecasting and Control, San
Francisco: Holden-Day, 1976.
Brenner, A. C., C. J. Koblinsky, and B. D. Beckley, A preliminary estimate of geoid-
induced variations in repeat orbit satellite altimeter observation, J. Geophys Res. 95,
3033–3040, 1990.
Brown, E. W., Theory of the motion of the Moon, Mem. Roy. Astron. Soc. 57, 51–145,
1905.
Calman, Y. Introduction to sea-surface topography from satellite altimetry, John Hopkins
APL Technical Digest, Laurel Md, 8(2), 206–211, 1987.
Carrère, L., and Lyard F., Modelling the barotropic response of the global ocean to
atmospheric wind and pressure forcing - comparisons with observations, Geophys. Res.
Lett., 30(6), 1275–1278, 2003.
Cartwright, D. E. and R. J. Taylor, New computations of tide-generating potential,
Geophys. J. Roy. Astron. Soc., 23, 45–74, 1971.
Cartwright, D. E. and A. C. Edden, Corrected Tables of Tidal Harmonics, Geophys. J.
Roy. Astron. Soc., 33, 253–264, 1973.
Cartwright, D. E. and R. D. Ray, Observations of the Mf ocean tide from Geosat
altimetry, Geophys. Res. Lett., 17(5), 619–622, 1990a.
Cartwright, D. E. and R. D. Ray, Oceanic Tides From Geosat Altimetry, J. Geophys.
Res., 95(C3), 3069–3090, 1990b.
Cartwright, D. E. and R. D. Ray, Energetics of global ocean tides from Geosat altimetry,
J. Geophys. Res., 96(C9), 16,897–16,912, 1991.
Cartwright, D. E. and R. D. Ray, On the radiational anomaly in the global ocean tide with
reference to satellite altimetry, Oceanologica Acta, 17(5), 453–459, 1994.

100
Casotto, S., Nominal ocean tide models for TOPEX/POSEIDON orbit determination,
Ph.D. dissertation, Center for Space Research, The University of Texas, Austin, 1989.
Chelton, D. B., J. C. Ries, B. J. Haines, L.-L. Fu, and P. Callahan, Satellite Altimetry, pp.
1–122. In: Satellite altimetry and earth sciences: A handbook of techniques and
applications, Lee-lueng Fu and Anny Cazenave (Eds.), Academic Press, 2001.
Chen, Y.Q., Recovery of terrestrial water storage change from low-low satellite-to-
satellite tracking, Rep. 485, Geodetic Science and Surveying, School of Earth Sciences,
Ohio State Univ., Columbus, OH, USA, 2007.
Cheng, C. C., Investigations into Green’s Function as Inversion-Free Solution of the
Kriging Equation, with Geodetic Applications, Rep. 472, Dept. of Civil and
Environmental Engineering and Geodetic Science, Ohio State Univ., Columbus, OH,
USA, 2004.
Cheng, Y. and O. B. Andersen, Multimission empirical ocean tide modeling for shallow
waters and polar seas, J. Geophys. Res., 116, C11001, doi:10.1029/2011JC007172, 2011.
Christensen, E. J., B. J. Haines, S. J. Keihm, C. S. Morris, R. A. Norman, G. H. Purcell,
B. G. Williams, B. D. Wilson, G. H. Born, M. E. Parke, S. K. Gill, C. K. Shum, B. D.
Tapley, R. Kolenkiewicz, and R. S. Nerem, Calibration of TOPEX/POSEIDON at
Platform Harvest, J. Geophys. Res., 99(C12), 24,465–24,485, 1994.
Cochrane, D., and G. Orcutt, Application of least squares regression to relationships
containing autocorrelated error terms, J. Am. Stat. Assoc., 44, 32–61, 1949.
Darwin, G. H., Report of a committee for the harmonic analysis of tidal observations,
Brit. Ass. Rep., 48–118, 1883.
Davis, R. W., J. G. Ortega-Ortiz, C. A. Ribic, W. E. Evans, D. C. Biggs, P. H. Ressler, R.
B. Cady, R. R. Leben, K. D. Mullin, and B. Würsig, Cetacean habitat in the northern
oceanic Gulf of Mexico. Deep-Sea Res. I 49: 121–142, 2001.
De Cesare, L. D.E Myers, and D. Posa, Estimating and modeling space-time correlation
structures, Statistics & Probability Letters, 51, 9–14, 2001.
De Cesare, L. D.E Myers, and D. Posa, FORTRAN programs for space-time modeling,
Computers and Geosciences, 28, 205–212, 2002.
Desai, S. D. and J. M. Wahr, Empirical ocean tide models estimated from
TOPEX/POSEIDON altimetry, J. Geophys. Res., 100(C12), 25,205–25,228, 1995.
Desai, S. D. Ocean tides from TOPEX/POSEIDON altimetry with some geophysical
applications, PhD Dissertation, Dept. of Aerospace Engineering Sciences, University of
Colorado, 1996.
Desai, S. D., Observing the pole tide with satellite altimetry, J. Geophys. Res., 107(C11),
3186, doi:10.1029/2001JC001224, 2002.
Desai, S. D. and B. J. Haines, Precise near-real-time sea surface height measurements
from the Jason-1 and Jason-2/OSTM missions, Mar. Geod., 33(S1): 419–434, 2010.
101
Ding X., and R. Coleman, Multiple outlier detection by evaluating redundancy
contributions of observations. J. Geod., 70, 489–498, 1996.
Doodson, A. T., The harmonic development of the tide-generating potential, Proceedings
of the Royal Society, 100, 305–329, 1921.
Dorandeu, J. and P. -Y. Le Traon, Effects of global mean atmospheric pressure variations
on mean sea level changes from TOPEX/POSEIDEN, Journal of Atmospheric and
Oceanic Technology, 16(9), 1279–1283, 1999.
Dorandeu, J., M. Ablain, and P. -Y. Le Traon, Reducing cross-track geoid gradient errors
around TOPEX/Poseidon and Jason-1 nominal tracks: Application to calculation of sea
level anomalies, Journal of Atmospheric and Oceanic Technology, 20(12), 1826–1838,
2003.
Eanes, R. and S. Bettadpur, The CSR3.0 global ocean tide model, Center for Space
Research, Technical Memorandum, CSR-TM-95-06, Univ. of Texas at Austin, 1995.
Egbert, G. D., A. F. Bennett, and M. G. G. Foreman, TOPEX/POSEIDON tides estimated
using a global inverse model, J. Geophys. Res., 99(C12), 24,821–24,852, 1994.
Egbert, G. D. and S. Y. Erofeeva, Efficient inverse modeling of the barotropic ocean
tides, Journal of Atmospheric and Oceanic Technology, 19, 183–204, 2002.
Escudier, P. and J. –L. Fellous. 2009. The Next 15 years of Satellite Altimetry: Ocean
surface topography constellation user requirements document. By Collecte Localization
Satellites (CLS) and Orange Bleue Conseil avaiable at
http://www.aviso.oceanobs.com/es/actualitades/news-detail/index.html?tx_ttnews
[tt_news]=609&tx_ttnews[backPid]=276&cHash=1a212e08df
Fang G., Y. Wang, Z. Wei, B. H. Choi, X. Wang, and J. Wang, Empirical cotidal charts
of the Bohai, Yellow, and East China Seas from 10 years of TOPEX/Poseidon altimetry.
J. Geophys. Res. 109: C11006, doi: 10.1029/2004JC002484, 2004.
Farrell, W. E., Deformation of the Earth by surface loads, Rev. Geophys. And Space
Phys., 10, 761–797, 1972.
Faugere Y., A. Ollivier, and J. F. Legeais, Envisat GDR quality assessement report
(cyclic), Cycle 095, Technical note SALP-RP-P2-EX-21121-CLS095, 2011 available at
http://www.aviso.oceanobs.com/html/donnees/calval/validation_report/en/welcome
_uk.hml
Fok, H. S., H. B. Iz, C. K. Shum, Y. Yi, O. Andersen, A. Braun, Y. Chao, G. Han, C. Y.
Kuo, K. Matsumoto, and Y. T. Song, Evaluation of ocean tide models used for Jason-2
altimetry corrections, Mar. Geod., 33(S1), 285–303, 2010.
Fotheringham, A. S., C. Brunsdon, and M. Charlton, Geographically weighted regression:
The analysis of spatially varying relationships, Wiley, Chichester, Hoboken, NJ, 2002.

102
Fotopoulos, G. , An analysis on the optimal combination of geoid, orthometric and
ellipsoidal height data, PhD thesis, Dept of Geomatics Engineering, University of
Calgary, Alberta, Canada, 2003.
Fu, L. L., E. J. Christensen, C. A. Yamarone Jr., M. Lefebvre, Y. M´enard, M. Dorrer,
and P. Escudier, TOPEX/POSEIDON mission overview, J. Geophys. Res., 99(C12),
24,369–24,381, 1994.
Fu L-L. and A. Cazenave, Satellite altimetry and earth sciences: a handbook of
techniques and applications. San Diego, Calif.; London: Academic Press, 2001.
Grafarend, E., A. Kleusberg, and B. Schaffrin, An introduction to the variance-
covariance-component estimation of Helmert type, Zeitschrift für Vermessungswesen,
105, 161–180, 1980.
Gross, R. S., The effect of ocean tides on the Earth's rotation as predicted by the results of
an ocean tide model, Geophys. Res. Lett., 20, 293-296, 1993.
Groves, G. W. and R. W. Reynolds, An orthogonalized convolution method of tide
prediction, J. Geophys. Res., 80(30), 4131–4138, 1975.
Guman, M. D., Determination of global mean sea level variations using multi-satellite
altimetry, Report CS-97-3, Center for Space Research, The University of Texas at Austin,
Austin, Texas, 1997.
Guo, J., J. Ou, and H. Wang, Robust estimation for correlated observations: two local
sensitivity-based downweighting strategies, J. Geod., 84, 243–250, 2010.
Guo, J. Y., O. Dierks, J. Neumeyer, and C. K. Shum, Weighting algorithms to stack
superconducting gravimeters data for the potential detection of the Slichter modes, J.
Geodyn., 41, 326–333, 2006.
Haines, B. J., G. H. Born, G. W. Rosborough, J. G. Marsh, and R. G. Williamson, Precise
orbit computation for the Geosat Exact Repeat Mission, J. Geophys. Res., 95(C3), 2871–
2885, 1990.
Han G., Three-dimensional modeling of tidal currents and mixing quantities over the
Newfoundland Shelf. J. Geophys. Res., 105(C5), 11407-11422, 2000.
Han, G., S. Paturi, B. deYoung, Y. Yi, and C.K. Shum, A 3-D data-assimilative tide
model of Northwest Atlantic. Atmosphere-Ocean, 48(1), 39–57, 2010.
Hampel F. R., E. M. Ronchetti, P. J. Rousseeuw, W. A. Stabel, Robust Statistics. The
approach based on influence functions. Wiley, New York. 1986.
Hekimoglu, S., and R. C. Erenoglu, Effect of heteroscedasticity and heterogeneousness
on outlier detection for geodetic networks, J Geod 81, 137–148, 2007.
Helmert, F. R., Die Ausgleichungsrechnung nach der Methode der kleinsten Quadrate,
Zweite Auflage, Teubner, Leipzig, Germany, 1907.

103
Helmert, F. R., Die Ausgleichsrechnung nach der Methode der kleinsten Quadrate, 3.
Aufl., Teubner, Leipzig, Germany, 1924.
Hendershott, M. C., Long Waves and Ocean Tides, In: Evolution in Physical
Oceanography, edited by B.A. Warren and C. Wunsch, MIT Press, pp. 292–341, 1981.
Hobson, E. W., The Theory of Spherical and Ellipsoidal Harmonics, Cambridge
University Press, 1965.
Holloway, G. and T. Sou, Has Arctic sea ice rapidly thinned? J. Climate, 15, 1691–1701,
2002.
Huber, P.J., Robust estimation of a location parameter, Annals Math. Stat. 35(1): 73–101,
1964.
Huber, P.J., Robust Statistics, Wiley, New York, 1981.
Huess, V. and O. B. Andersen, Seasonal variation in the main tidal constituent from
altimetry. Geophys. Res. Lett., 28(4) 567–570, 2001.
IERS Conventions, Gérard Petit and Brian Luzum (eds.). (IERS Technical Note ; 36)
Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie, 179 pp.,
2010.
Inazu, D., T. Sato, S. Miura, Y. Ohta, K. Nakamura, H. Fujimoto, C. F. Larsen, and T.
Higuchi, Accurate ocean tide modeling in southeast Alaska and large tidal dissipation
around Glacier Bay. J. Oceanogr. 65: 335–347, 2009.
Iz, H. B., and Y. Q. Chen, VLBI rates with first order autoregressive disturbances, J.
Geodyn., 28, 131–145, 1999.
Iz, H. B., Polar motion modeling, analysis, and prediction with time dependent harmonic
coefficients, J. Geod., 82, 871–881, 2008.
Ito, T., and M. Simons, Probing asthenospheric density, temperature, and elastic moduli
below the western United States, Science, 332, 947–951, 2011.
Johnson R. A., and D. W. Wichern, Applied Multivariate Statistical Analysis, 3rd edn,
Prentice-Hall, Englewood Cliffs, New Jersey, 1992.
Jordan, S. K., Self-consistent statistical models for the gravity anomaly, vertical
deflections, and undulation of the geoid, J. Geophys. Res, 77, 3660–3670, 1972.
Jørgensen, P. C., Kubik, K., Frederiksen, P., and W. Weng, Ah, robust estimation!.
Australian Journal of Geodesy, Photogrammetry, and Surveying, 42, 19–32, 1985.
Junkins, J. L., G. W. Miller, and J. R. Jancaitis, A weighting function approach to
modeling of irregular surfaces, J. Geophys. Res, 78, 1794–1803, 1973.
Kagan, B. A., D. A. Romanenkov, and E. V. Sofina, Tidal ice drift and ice-generated
changes in the tidal dynamics/energetic on the Siberian continental shelf. Continental
Shelf Research, 28, 351–368, 2008.

104
Kagan, B. A., E. V. Sofina, and A. A. Timofeev, Modeling of the M2 surface and internal
tides and their seasonal variability in the Arctic Ocean: Dynamics, energetic and tidally
induced diapycnal diffusion. Journal of Marine Research, 69, 245–276, 2011.
Kaula, W. M., The terrestrial environment: solid earth and ocean physics, NASA Rep.
study at Williamstown, MA, NASA CR-1579, 1969.
Kavouras, M., On the detection of outliers and the determination of reliability in geodetic
network, Technical Report No. 87, Dept. of Geodesy and Geomatics Engineering,
University of New Brunswick, Fredericton, N. B., Canada, 1982.
Kern, M., T. Preimesberger, M. Allesch, R. Pail, J. Bouman, and R. Koop, Outlier
detection algorithms and their performance in GOCE gravity field processing, J. Geod.,
78, 509–519, 2005.
King, C., D. Stammer, and C. Wunsch, Tide model comparison at CPMO/MIT. Working
paper for the TOPEX/POSEIDON science working team tide model study. Dept. of
Earth, Atmos., and Planet. Sc., MIT, Massachusetts, 1995.
King, M., N. T. Penna, P. J. Clarke, and E. C. King, Validation of tide models around
Antarctica using onshore GPS and gravity data, J. Geophys. Res., 110, B08401,
doi:10.1029/2004JB003390, 2005.
Kleusberg, A., and E. W. Grafarend, Expectation and variance component estimation of
multivariate gyrotheodolite observation II, Allg. Vermessungsnachrichten 88, 104–108,
1981.
Knight, N. L. and J. Wang, A comparison of outlier detection procedures and robust
estimation methods in GPS positioning, J. Navigation, 62, 699–709, 2009.
Knudsen, P., Global low harmonic degree models of the seasonal variability and residual
ocean tides from TOPEX/POSEIDON altimeter data, J. Geophys. Res., 99(C12), 24,643–
24,655, 1994.
Koch, K. R., Outlier tests and reliability measures. Vermessungswesen und
Raumordnung, 45, 400–411, 1983 (in German).
Koch, K. R., Parameter estimation and hypothesis testing in linear models, 2nd edition,
Springer Berlin, 1999.
Kotsakis, C., and M. G. Sideris, On the adjustment of combined GPS/leveling/geoid
networks, J. Geod., 73, 412–421, 1999.
Krarup, T., K. Kubik, and J. Juhl, Gӧtterdaemmerung over least squares adjustment, In:
Proceedings of International Society for Photogrammetry 14th Congress, Vol. 33,
Hamburg, 369–378, 1980.
Kruizinga, G., Validation and applications of satellite radar altimetry. PhD dissertation.
Center of Space Research, University of Texas at Austin, 1997.
Kubik, K., W. Wang, and P. Frederiksen, Oh, gross errors!. Australian Journal of
Geodesy, Photogrammetry, and Surveying, 42, 1–18, 1985.
105
Kudryavtsev, S. M., Improved harmonic development of the Earth tide-generating
potential. J. Geod., 77(12), 829–838, 2004.
Kusche, J., A Monte-Carlo technique for weight estimation in satellite geodesy, J. Geod.,
76(11–12), 641–652, 2003.
Kwok, R., and N. Untersteiner, The thinning of Arctic sea ice. Phys. Today, 64(4), 36-41,
2011.
Kwon, J. H., D. Grejner-Brzezinska, T.-S. Bae, and C.-K. Hong, A triple difference
approach to low Earth orbiter precision orbit determination, Journal of Navigation, 56
(3), 457–473, 2003.
Kyriakidis, P. C., A. G. Journel, Geostatistical space-time models: A review, Math Geol.,
31 (6), 651–684, 1999.
Lambeck, K., Geophysical Geodesy, The Slow Deformation of the Earth, Oxford
University Press, 1988.
Laplace, P. S., Recherches sur quelque points du systeme du monde, Mem. Acad. Roy.
Soc., 75–182, 1775.
Latychev K., J. X. Mitrovica, M. Ishii, N. Chan, and J.L. Davis, Body tides on a 3–D
elastic earth: Toward a tidal tomography, Earth Planet. Sci. Lett., 277, 86–90, 2009.
Leben, R. R., G. Born, J. D. Thompson, and C. A. Fox, Mean sea surface and variability
of the Gulf of Mexico using geosat altimetry data. J. Geophys. Res. 95 (C3): 3025–3032,
1990.
Lee, H. K., Radar altimetry methods for solid earth geodynamics studies, Rep. 489,
Geodetic Science and Surveying, School of Earth Sciences, Ohio State Univ., Columbus,
OH, USA, 2008.
Leeuwenburgh, O., O. B. Andersen, and V. Huess, Seasonal tide variations from tide
gauges and altimetry. Phys. Chem. Earth (A), 24(4), 403-406, 1999.
Lefevre, F. F. H. Lyard, and C. Le Provost, FES98: A new global tide finite element
solution independent of altimetry, Geophys. Res. Lett., 27(17), 2717–2720, 2000.
Lemoine, F. G., N. P. Zelensky, D. S. Chinn, B. D. Beckley, and J. L. Lillibridge,
Towards the GEOSAT Follow-On precise orbit determination goals of high accuracy and
near-real-time processing, AIAA/ASS Astrodynamics Specialist Conference, Keystone,
Colorado, August 21–24, 2006.
Le Provost, C. and P. Vincent, Extensive test of precision for a finite element model of
ocean tide, J. Comput. Phys., 65, 273–291, 1986.
Le Provost, C., M. L. Genco, F. Lyard, P. Vincent, and P. Canceil, Spectroscopy of the
world ocean tides from a finite element hydrodynamic model, J. Geophys. Res., 99(C12),
24,777–24,797, 1994.

106
Le Provost, C., Ocean tides, In: Satellite altimetry and earth sciences, Lee-lueng Fu and
Anny Cazenave (Eds.), International Geophysics Series, Vol. 69, Academic Press, 267–
303, 2001.
Le Traon, P. Y., F. Nadal, and N. Ducet, An improved mapping method of multisatellite
altimeter data, J. Atmos. Oceanic Technol., 15, 522–534, 1998.
Lorell, J., E. Colquitt, and R. J. Anderle, Ionospheric correction for SEASAT altimeter
height measurement, J. Geophys. Res., 87(C5), 3207–3212, 1982.
Lyard F., F. Lefevre, T. Letellier, and O. Francis, Modelling the global ocean tides:
Modern insights from FES2004. Ocean Dyn. 56 (5–6): 394–415, 2006.
Matsumoto, K., T. Takanezawa, and M. Ooe, Ocean tide models developed by
assimilating TOPEX/POSEIDON altimeter data into hydrodynamic model: A global
model and a regional model around Japan, J. Oceangr., 56, 567–581, 2000.
Menard, Y., L.-L. Fu, P. Escudier, F. Parisot, J. Perbos, P. Vincent, S. Desai, B. Haines,
and G. Kunstmann, The Jason-1 mission, Mar. Geod., 26(3), 131–146, 2003.
Moritz, H., Least-squares collocation, Rev. Geophys. Space Phys., 16, 421–430, 1978.
Munk, W. H. and D. E. Cartwright, Tidal spectroscopy and prediction, Phil. Trans. Roy.
Soc. A, 259(1105), 553–581, 1966.
Munk, W. H. and G. J. F. MacDonald, The Rotation of the Earth, Cambridge University
Press, 1975.
Newcomb, S., Astronomical constants (The elements of the four inner planets and the
fundamental constants of astronomy), Supplement to the American Ephemeris and
Nautical Almanac for 1897, US Government Printing Office, Washington, DC, 1895.
Niedzielski, T. and W. Kosek, Forcasting sea level anomalies from TOPEX/Poseidon and
Jason-1 satellite altimetry, J. Geod., 83, 469–476, 2009.
Padman, L. and S. Erofeeva, A Barotropic Inverse Tidal Model for the Arctic Ocean,
Geophys. Res. Lett., 31(2), L02303, doi:10.1029/2003GL019003, 2004.
Parke, M. E., R.H. Stewart, D. L. Farless, and D. E. Cartwright, On the choice of orbits
for an altimetric satellite to study ocean circulation and tides, J. Geophys. Res., 92(C11),
11,693–11,707, 1987.
Pascual, A., M. Marcos, and D. Gomis, Comparing the sea level response to pressure and
wind forcing of two barotropic models: Validation with tide gauge and altimetry data, J.
Geophys. Res., 113, C07011, doi:10.1029/2007JC004459, 2008.
Ponte, R. M., D. A. Salstein and R. D. Rosen, Sea level response to pressure forcing in a
barotropic numerical model, J. Phys. Oceanog., 21, 1043-1057, 1991.
Pope, A., The statistics of residuals and detection of outliers, NOAA Technical Rep. No.
66, National Geodetic Survey, Rockville, Maryland, 1976.
Pugh, D. T., Tides, Surges and Mean Sea-Level, John Wiley & Sons, 1987.
107
Pugh, D., Changing Sea Levels: Effects of Tides, Weather and Climate, Cambridge
University Press, 2004.
Rangelova, E., G. Fotopoulos, and M. G. Sideris, On the use of iterative reweighting
least-squares and outlier detection for empirically modeling rates of vertical
displacement, J. Geod., 83, 523–535, 2009.
Ray, R. D., Global ocean tide models on the eve of TOPEX/POSEIDON, IEEE Trans.
Geosci. Remote Sens., 31(2), 355–364, 1993.
Ray, R. D., D. J. Steinberg, B. F. Chao, and D. E. Cartwright, Diurnal and semidiurnal
variations in the Earth's rotation rate induced by oceanic tides, Science, 264, 830–832,
1994.
Ray, R.D., R. J. Eanes, and B. F. Chao, Detection of tidal dissipation in the solid earth by
satellite tracking and altimetry, Nature, 381, 595–597, 1996.
Ray, R.D., Ocean self-attraction and loading in numerical tidal models, Mar. Geod., 21,
181–192, 1998.
Ray, R. D., A global ocean tide model from TOPEX/POSEIDON altimetry: GOT99.2.
NASA Tech Memo – 209478, 1999.
Ray, R. D., G. D. Egbert, and S. Y. Erofeeva, Tidal predictions in shelf and coastal
waters: Status and prospects. In: Coastal Altimetry, Springer, 191–216, 2011.
Renganathan, V., Arctic sea ice freeboard heights from satellite altimetry. PhD
dissertation, Dept. of Geomatics Engineering, University of Calgary, 2010.
Rodriguez-Iturbe, I., and J. M. Meija, The design of rainfall networks in time and space.
Water Resources Research 10, 713–728, 1974.
Rouhani, S., and T. J. Hall, Space-time kriging of ground water data. Proceedings 3rd
International Geostatistics Congress. (Avignon, France), Kluwer Academic Publishers,
Dordrecht, Vol. 2, 639–651, 1989.
Rousseeuw, P. J., Leroy A.M., Robust Regression and Outlier Detection. Wiley, New
York, 1987.
Rudenko, S., M. Otten, P. Visser, R. Scharroo, T. Schone, and S. Esselborn, New
improved orbit solutions for the ERS-1 and ERS-2 satellites, Adv. Space Res., 49, 1229–
1244, 2012.
Rummel, R., Principle of satellite altimetry and elimination of radial orbit errors, In:
Satellite altimetry in geodesy and oceanography, Lecture Notes in Earth Sciences, 50,
191–241, 1993.
Saastamoinen, J., Atmospheric correction for the troposphere and stratosphere in radio
ranging of satellites, Geophysical Monograph 15, American Geophysical Union, 1972.
Sabaka, T. J., D. D. Rowlands, S. B. Luthcke, and J. -P. Boy, Improving global mass flux
solutions from Gravity Recovery and Climate Experiment (GRACE) through forward

108
modeling and continuous time correlation, J. Geophys. Res., 115, B11403,
doi:10.1029/2010JB007533, 2010.
Sanchez, B. V. and N. L. Pavlis, Estimation of main tidal constituents from TOPEX
altimetry using a Proudman function expansion, J. Geophys. Res., 100(C12), 25,229–
25,248, 1995.
Savcenko, R. and W. Bosch, EOT08a – empirical ocean tide model from multi-mission
satellite altimetry. In: Report No. 81 Deutsches Geodätisches Forschungsinstitut (DGFI),
München, Germany, 2008.
Savcenko, R. and W. Bosch, EOT10a – a new result of empirical ocean tide modeling,
Ocean Surface Topography Science Team (OSTST) Meeting, 18-20 October, Lisbon,
Portugal, 2010.
Savcenko, R. and W. Bosch, EOT11a – a new tide model from multi-mission satellite
altimetry, Ocean Surface Topography Science Team (OSTST) Meeting, 19-21 October,
San Diego, USA, 2011.
Schaffrin, B., Varianz-Kovarianz-Komponenten-Schätzung bei der Ausgleichung
heterogener Wiederholungsmessungen, Publ. DGK C-282, Deutsche Geodätische
Kommission, München, 1983.
Schaffrin, B., Reliability measures for correlated observations, J Surv Eng ASCE, 123(3),
126–137, 1997.
Schaffrin, B., and H. B. Iz, Integrating heterogeneous data sets with partial
inconsistencies, in: M.G.Sideris (ed.), Gravity, Geoid and Geodynamics 2000, Springer
Series, IAG-Symp., Berlin, 123, 49–54, 2001.
Schrama, E. J. O. and R. D. Ray, A preliminary tidal analysis of TOPEX/POSEIDON
altimetry, J. Geophys. Res., 99(C12), 24,799–24,808, 1994.
Schureman, P., Manual of Harmonic Analysis and Prediction of Tides, U.S. Coast and
Geodetic Survey, Special Publication No. 98, 1971.
Schwiderski, E. W., Ocean Tides, Part I: Global Ocean Tidal Equations, Mar. Geod., 3,
161–217, 1980.
Schwiderski, E. W., Atlas of Ocean Tidal Charts and Maps, Part I: The Semidiurnal
Principal Lunar Tide M2, Mar. Geod., 6(3–4), 219–265, 1983.
Seeber, G., Satellite Geodesy, Walter de Gruyter, Berlin, 1993.
Sheinbaum J., J. Candela, A. Badan, and J. Ochoa, Flow structure and transport in the
Yucatan Channel. Geophys. Res. Lett. 29 (3): 10.1029/2001GL013990, 2002.
Shum, C. K., P. L. Woodworth, O. B. Andersen, G. D. Egbert, O. Francis, C. King, S. M.
Klosko, C. Le Provost, X. Li, J. M. Molines, M. E. Parke, R. D. Ray, M. G. Schlax, D.
Stammer, C. C. Tierney, P. Vincent, and C. I. Wunsch, Accuracy assessment of recent
ocean tide models, J. Geophys. Res., 102(C11), 25,173–25,194, 1997.

109
Shum, C. K., N. Yu, and C. S. Morris, Recent advances in ocean tidal science. J. Geod.
Soc. Japan, 47(1), 528–537, 2001.
Shum, C. K., Y. C. Yi, H. K. Li, K. Matsumoto, T. Sato, X. C. Wang, Y. Chao, X. L.
Deng, H. B. Iz, Coastal Ocean Tide Modeling Using Satellite Altimetry. Ocean Surface
Topography Science Team (OSTST) Meeting, Venice, Italy, 2006.
Sisman, Y., Outlier measurement analysis with the robust estimation, Scientific Research
and Essay, 5(7), 668–678, 2010.
Smith, A. J. E., Application of satellite altimetry for global ocean tide modeling. PhD
Dissertation. Delft Institute for Earth-Oriented Space Research, Delft University Press,
1999.
Stewart, R. H., Introduction to Physical Oceanography, 2007,
http://oceanworld.tamu.edu/resources/ocng_textbook/contents.html.
Taff, L. G. , Celestial Mechanics, A Computational Guide for the Practitioner, John
Wiley & Sons, 1985.
Tapley, B. D., J. B. Lundverg, and G. H. Born, The SEASAT altimeter wet tropospheric
range correction, J. Geophys. Res., 87, 3213–3220, 1982.
Tapley, B. D., and G. W. Rosborough, Geographically correlated orbit error and its effect
on satellite altimetry, J. Geophys. Res., 90(C6), 11,817–11,831, 1985.
Tapley, B. D., J. C. Ries, G. W. Davis, R. J. Eanes, B. E. Schutz, C. K. Shum, M. M.
Watkins, J. A. Marshall, R. S. Nerem, B. H. Putney, S. M. Klosko, S. B. Luthcke, D.
Pavlis, R. G. Williamson, and N. P. Zelensky, Precise orbit determination for
TOPEX/POSEIDON, J. Geophys. Res., 99(C12), 24,383–24,404, 1994a.
Tapley, B. D., D. P. Chambers, C. K. Shum, R. J. Eanes, J. C. Ries, and R. H. Stewart,
Accuracy assessment of the large-scale dynamic ocean topography from
TOPEX/POSEIDON altimetry, J. Geophys. Res., 99(C12), 24,605–24,617, 1994b.
Van Loon, J.P., Functional and stochastic modeling of satellite gravity data. PhD
Dissertation. Delft Institute for Earth-Oriented Space Research, Delft University Press,
2008.
Wagner, C. A., Summer school lectures on satellite altimetry, Lecture Notes in Earth
Sciences, Theory of satellite geodesy and gravity field determination, 25, Springer-
Verlag, 1989.
Wahr, J. M., Body tides on an elliptical, rotating, elastic, oceanless Earth, Geophys. J. R.
Astr. Soc., 64, 677–703, 1981.
Wahr, J. M. and T. Sasao, A diurnal resonance in the ocean tide and in the Earth’s load
response due the resonant free “core nutation”, Geophys. J. R. Astr. Soc., 64, 747–765,
1981.
Wahr, J.M., Deformation induced by polar motion, J. Geophys. Res., 90, 9363–9368,
1985.
110
Wang, Y., Ocean tide modeling in the Southern Ocean, Rep. 471, Geodetic Science and
Surveying, School of Earth Sciences, Ohio State Univ., Columbus, OH, USA, 2004.
Wenzel, H.- G., Tide-generating potential for the Earth, Lecture Notes in Earth Sciences,
66, 9–26, 1997.
Yang, Y., Robust estimation of geodetic datum transformation, J. Geod., 73, 268–274,
1999.
Yang, Y., Xu, T., and Song L., Robust estimation of variance components with
application in global positioning system network adjustment. J. Surv. Eng. ACSE, 131(4),
107–112, 2005.
Yi, Y., Determination of gridded mean sea surface from TOPEX, ERS-1 and GEOSAT
altimeter data, Rep. 434, Dept. Geod. Sci. Surv., The Ohio State Univ., Columbus, Ohio,
1995.
Yi, Y., K. Matsumoto, C. K. Shum, Y. Wang, and R. Mautz, Advances in southern ocean
tide modeling. J. Geodyn., 41, 128–132, 2006.
Yi, Y., The Ohio State University stackfiles for satellite radar altimeter data, Rep. 495,
Geodetic Science and Surveying, School of Earth Sciences, Ohio State Univ., Columbus,
Ohio, USA, 2010.
Yuan, L. G., X. L. Ding, P. Zhong, W. Chen, and D. F. Huang. 2009. Estimates of ocean
tide loading displacements and its impact on position time series in Hong Kong using a
dense continuous GPS network. J. Geod. doi: 10.1007/s00190-009-0319-0.
Zu, T., J. Gan, and S. Y. Erofeeva, Numerical study of the tide and tidal dynamics in the
South China Sea. Deep-Sea Research I, 55, 137–154, 2008.

111
 Appendix A: Tables for Overall Results

Table A.1 Ground truth comparison of ocean tide models at 102 pelagic tide gauge stations (in cm). Note
that D is the percentage value representing the discrepancy of the ocean tide models against the ground
truth data (i.e.RSS/RSSIQ×100%). Bold values indicate the models with the best performance in this test.
Note that bold values indicate the best model solution over a particular tidal constituent and/or the overall
discrepancy in terms of RSS and D (%).

Model M2 S2 K1 O1 N2 P1 K2 Q1 RSS D (%)

RSSIQ = 39.01 cm
DTU10 1.39 0.92 0.98 0.74 1.52 0.82 1.09 0.39 2.93 7.51
EOT11a 1.43 0.94 0.95 0.73 1.54 0.80 1.08 0.38 2.95 7.56
EOT10a 1.41 0.93 0.96 0.73 1.55 0.80 1.07 0.37 2.94 7.53
EOT08a 1.44 1.06 0.98 0.73 1.57 0.84 1.11 0.40 3.04 7.78
FES2004 1.45 0.95 1.00 0.75 1.56 0.83 1.10 0.40 3.01 7.71
GOT00.2 1.45 1.03 1.01 0.85 1.51 0.77 1.05 0.39 3.00 7.70
GOT4.7 1.46 1.01 1.01 0.76 1.54 0.76 1.05 0.37 2.99 7.67
TPXO.6.2 1.46 0.92 1.04 0.81 1.56 0.70 1.06 0.36 2.98 7.64
TPXO.7.1 1.41 0.92 1.07 0.79 1.57 0.70 1.08 0.37 2.98 7.61
TPXO.7.2 1.43 0.92 1.07 0.86 1.57 0.69 1.07 0.37 3.00 7.66
OSU12sw 1.30 1.12 0.91 0.72 1.57 0.71 1.11 0.44 2.94 7.55
OSU12vce 1.30 1.15 0.92 0.73 1.56 0.71 1.10 0.44 2.95 7.56

Table A.2 Ground truth comparison of ocean tide models with common tide gauges at 739 coastal tide
gauge stations (in cm).

M2 S2 K1 O1 N2 P1 K2 Q1
Model RSS D (%)
(548) (549) (546) (541) (487) (501) (505) (367)
RSSIQ = 75.06 cm
DTU10 18.27 8.51 7.43 4.57 5.31 2.62 2.89 1.25 22.96 30.58
EOT11a 20.38 9.23 8.52 4.61 5.38 3.16 3.10 1.24 25.38 33.82
EOT10a 20.61 9.37 8.61 4.98 5.51 3.12 3.11 1.26 25.75 34.30
EOT08a 20.06 9.10 8.16 4.78 5.36 3.04 3.08 1.20 24.97 33.26
FES2004 21.73 9.73 8.86 5.22 5.64 3.16 3.32 1.29 26.96 35.92
GOT00.2 17.16 8.52 6.10 3.98 4.89 2.27 3.05 1.16 21.44 28.57
GOT4.7 17.08 8.34 6.11 3.74 4.80 2.21 3.12 1.11 21.25 28.31
TPXO.6.2 25.20 11.37 5.53 4.11 7.18 2.12 3.67 1.39 29.72 39.59
TPXO.7.1 25.38 14.61 5.57 3.84 6.84 2.80 4.81 3.89 31.51 41.98
TPXO.7.2 24.82 17.61 5.64 4.94 7.06 2.29 4.21 2.00 32.54 43.35
OSU12sw 15.95 7.85 6.20 3.87 4.98 2.25 2.62 1.62 20.21 26.92
OSU12vce 16.88 8.28 5.74 3.79 4.95 2.07 2.73 1.57 20.96 27.92

112
 Table A.3 Standard deviations of residual SSH anomaly of ocean tide models along satellite tracks in the overall ocean. Stdev (after) are the standard
deviation of the SSH anomaly after ocean tide correction; RSS Stdev (before) and RSS Stdev (after) are the Root Sum of Squares of the standard deviation
of the SSH anomaly before and after ocean tide correction for the entire altimeter data, respectively; VE is the variance explained by ocean tide correction.
Note that bold values indicate the best model solution over a particular altimeter data and/or the overall variance reduction, and Jason-2 altimeter data were
not included in the generation of the OSU12 model solution.

Stdev (after) RSS

Model TOPEX/ TOPEX Jason-1 Stdev VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
34.70 34.55 34.88 33.19 34.15 33.87 35.24
(before)
DTU10 10.61 9.52 10.55 10.31 11.54 9.87 11.41 27.96 69.26
EOT11a 10.62 9.51 10.56 10.28 11.52 9.88 11.39 27.94 69.28
EOT10a 10.64 9.52 10.57 10.30 11.55 9.87 11.41 27.97 69.24
113

EOT08a 10.65 9.54 10.56 10.30 11.55 9.89 11.43 28.00 69.21
FES2004 10.71 9.61 10.63 10.37 11.60 9.95 11.51 28.17 69.03
GOT00.2 10.64 9.57 10.58 10.32 11.56 9.90 11.44 28.03 69.18
GOT4.7 10.62 9.52 10.54 10.29 11.54 9.88 11.41 27.95 69.26
TPXO.6.2 10.71 9.64 10.63 10.38 11.58 9.94 11.49 28.16 69.03
TPXO.7.1 10.67 9.59 10.60 10.37 11.55 9.91 11.46 28.09 69.11
TPXO.7.2 10.65 9.57 10.59 10.35 11.56 9.91 11.44 28.04 69.16
OSU12sw 10.80 9.33 10.40 10.20 11.64 9.77 11.52 27.91 69.31
OSU12vce 10.78 9.32 10.40 10.19 11.61 9.76 11.49 27.87 69.35
RSS Stdev (before) = 90.94 cm
 Table A.4 Standard deviations of residual SSH anomaly of ocean tide models along satellite tracks in the in Shallow ocean with depth less than 1000 m (in
cm).

Stdev (after) RSS

Model TOPEX/ TOPEX Jason-1 Stdev VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
54.14 53.61 54.23 52.46 53.08 53.51 54.68
(before)
DTU10 11.18 10.52 11.36 11.33 12.02 10.46 12.19 29.93 78.92
EOT11a 11.37 10.68 11.46 11.38 12.04 10.56 12.34 30.22 78.72
114

EOT10a 11.45 10.74 11.53 11.42 12.13 10.60 12.42 30.39 78.60
EOT08a 11.41 10.72 11.48 11.41 12.13 10.61 12.44 30.35 78.63
FES2004 12.07 11.29 12.10 11.86 12.56 11.36 13.08 31.91 77.53
GOT00.2 11.47 10.97 11.65 11.51 12.25 10.70 12.49 30.66 78.41
GOT4.7 11.20 10.57 11.27 11.25 12.06 10.49 12.26 29.94 78.92
TPXO.6.2 12.07 11.77 12.26 12.17 12.83 11.31 13.08 32.34 77.23
TPXO.7.1 11.88 11.41 12.04 12.15 12.45 11.02 12.93 31.74 77.65
TPXO.7.2 11.56 11.09 11.75 11.75 12.25 10.78 12.60 30.95 78.21
OSU12sw 11.42 10.03 10.93 11.02 12.20 10.21 12.43 29.66 79.12
OSU12vce 11.42 10.07 10.97 11.05 12.14 10.19 12.38 29.64 79.13
RSS Stdev (before) = 142.02 cm
 Table A.5 Standard deviations of residual SSH anomaly of ocean tide models along satellite tracks in the in Deep Ocean with depth greater than 1000 m (in
cm).

Stdev (after) RSS

Model TOPEX/ TOPEX Jason-1 Stdev VE (%)
GFO Envisat Jason-1 Jason-2 (after)
Poseidon Interleave Interleave
Stdev
32.78 32.72 32.96 31.22 33.02 32.38 33.07
(before)
DTU10 10.57 9.44 10.49 10.23 11.52 9.84 11.34 27.81 67.75
EOT11a 10.57 9.42 10.49 10.20 11.49 9.84 11.30 27.77 67.81
EOT10a 10.57 9.43 10.49 10.21 11.52 9.82 11.32 27.79 67.78
EOT08a 10.59 9.45 10.48 10.21 11.52 9.85 11.34 27.82 67.74
FES2004 10.60 9.47 10.51 10.24 11.55 9.87 11.37 27.88 67.67
GOT00.2 10.58 9.46 10.49 10.23 11.53 9.86 11.35 27.83 67.73
GOT4.7 10.57 9.44 10.48 10.21 11.52 9.84 11.34 27.80 67.76
115

TPXO.6.2 10.60 9.47 10.50 10.23 11.52 9.86 11.34 27.84 67.71
TPXO.7.1 10.58 9.45 10.49 10.22 11.51 9.85 11.33 27.81 67.76
TPXO.7.2 10.57 9.45 10.49 10.23 11.51 9.85 11.34 27.82 67.74
OSU12sw 10.75 9.27 10.36 10.13 11.61 9.75 11.43 27.78 67.79
OSU12vce 10.73 9.27 10.35 10.12 11.58 9.74 11.41 27.74 67.84
RSS Stdev (before) = 86.24 cm