CHAPTER 2 MARINE GRAVITY ANOMALIES FROM SATELLITE
2.5 Multi-Satellite Altimeter Data Processing
2.5.3 Choice of Ocean Tide Model
Ocean tide creates a deviation of the instantaneous sea surface from the mean sea surface. There are now more than 10 global ocean tide models (Shum et al., 1997) available for correcting tidal effect in altimetry. Table 2.4, partly from Matsumoto et al. (2000), shows the RMS collinear differences of T/P and the errors in along-track DOV using NAO99b (Matsumoto et al., 2000), CSR4.0 (Eanes, 1999) and GOT99.2b
Figure 2.3: Estimated standard deviations of point SSH of Geosat/ERM (top), ERS-1 (center) and TOPEX/POSEIDON (bottom)
(Matsumoto et al., 2001) tide models. Both NAO99b and GOT99.2b are based on hydrodynamic solutions and a further enhancement by assimilating T/P altimetry data into the solutions. The CSR tide models (3.0 and 4.0 versions) use the orthotide approach to model the residual tides of some preliminary hydrodynamic oceans tide
models using T/P altimeter data. From Table 2.4, it seems that NAO99b is the best
model among the three. Using T/P SSHs, Chen (2001) found that the NAO99b tide model yields the smallest RMS crossover differences of SSH compared to the CSR4.0 and GOT99.2b tide models. Over shallow waters, all the collinear differences exceed 10 cm, which translate to a 48-µrad (10-6 radian) error in DOV of Geosat/GM. Even in the deep oceans, the collinear difference-implied DOV errors are still very large.
Using available resources here, we conducted a test over the SCS to compare the accuracies of predicted gravity anomalies using the CSR3.0, CSR4.0 and GOT99.2b tide models. As shown in Table 2.5, the NAO99b tide model produces the best accuracy in gravity anomaly. However, as seen in Table 2.5, the differences in accuracy are very close. This should be partly due to the fact that DOV is insensitive to long wavelength tide model error. For Geosat/ERM, ERS-1/35d, ERS-2/35d and T/P, the use of CSR3.0, CRS4.0 and NAO99b will probably not make too much difference because of the reduction of tide model error by averaging data from repeat cycles.
Table 2.4: RMS collinear differences (in cm) of T/P (Matsumoto et al. 2000) and corresponding error (in µrad) in along-track DOV of Geosat/GM using different tide models
0 < Ha < 0.2 0.2 < H < 1 1 < H Tide model
Difference Error Difference Error Difference Error
NAO99b 11.20 45 6.98 28 8.56 35
CSR4.0 15.77 64 7.37 30 8.55 35
GOT99.2b 13.99 57 7.37 30 8.65 35
aH: ocean depth in km
Table 2.5: RMS differences (in mgals) between shipborne and altimeter-derived gravity anomalies with different tide models
Tide model ERS-1/GM Geosat/GM (JGM3)
NAO99b 11.90 9.77
CSR3.0 11.93 Not available
CSR4.0 Not available 9.78
CHAPTER 3
GLOBAL MODELS OF MEAN SEA SURFACE AND GRAVITY ANOMALY
3.1 Introduction
Satellite altimetry has begun a new era in earth sciences research. Mean sea surface height (MSSH) and marine gravity anomaly are the most two important products of satellite altimetry for geodetic and geophysical applications. MSSH is useful in numerous applications such as global tide modeling, sea level change study, reduction of altimeter observations to reference tracks, and bathymetry prediction.
Applications of marine gravity have been illustrated in, e.g., Sandwell and Smith (1997). Ignoring the sea surface topography (SST), marine gravity anomaly is basically equivalent to MSSH, and the two have a simple, linear relationship in the spectral domain. Current global MSSH models have been largely constructed by direct gridding of altimeter observed sea surface heights (SSHs), e.g., the MSSH models of Kort-og Matrikelstyrelsen (KMS) (Andersen and Knudsen, 1998), National Aeronautics and Space Administration (NASA)/Goddard Space Flight Center (GSFC) (Wang, 2000) and Collecte Localisation Satellites (CLS) (Hernandez and Schaeffer, 2000). A somewhat unconventional approach was employed by Yi (1995), who used a combined SSH and geoid gradient from multi-satellite missions to construct the OSU95 MSSH model, see also Rapp and Yi (1997). When using the direct gridding method, SSHs must be carefully crossover adjusted to remove inconsistency of SSHs at crossover points. Insufficient or improper crossover points
will easily lead to artifacts in the resulting field such as track pattern and extremely large signature. With MSSH computed, marine gravity anomaly can be obtained by a simple conversion in the spectral domain, see, e.g., Schwarz et al. (1990) for a complete derivation of the spectral conversion between MSSH and gravity anomaly (assuming that SST is removed).
As discussed in (2.2.3), using deflections of the vertical (DOV) as the data type is one of the choices in gravity anomaly recovery from satellite altimetry data. The major argument of using DOV is that DOV is less contaminated by long wavelength errors than SSH, and using DOV requires no crossover adjustment; see, e.g., Sandwell and Smith (1997), Hwang et al. (1998), and Andersen and Knudsen (1998). With DOV from altimetry, it is possible to compute MSSH using the deflection-geoid formula derived by Hwang (1998). Taking the advantage of DOV, in this chapter we will compute simultaneously a global MSSH grid and a global gravity anomaly grid by the deflection-geoid and the inverse Vening Meinesz formulae from multi-satellite altimetry data. In parallel to this study is the latest models of ocean tide and other geophysical corrections.
3.2 Forming north and east components of DOV
3.2.1 Computing along-track DOV
The methods for computing MSSH and gravity anomaly in this paper will use DOV as the data type. By definition, an along-track DOV is the geoid gradient with an opposite sign:
s g N
∂
−∂
= (3.1)
where N is the geoid, which is a surface function, and s is the along-track distance. By this definition one would need to first construct a surface of the geoid and then perform directional derivative along s to get DOV. Following this concept, we first fit a cubic spline (De Boor, 1978) to the along-track geoidal heights from altimetry. Then the along-track derivative is obtained by differentiating the spline. The actual numerical computations were done by the International Mathematical and Statistical Library (IMSL) routines. It turns that such a procedure results in very noisy DOV, which should be due to the interpolation error in fitting the spline. The interpolation error is particularly large when point intervals along a track segment are not uniform.
Although this approach seems rigorous, it did not produce good results. A better result is obtained by simply approximating DOV by the slope of two successive geoidal heights:
where d is the point spacing. The geographic location of g is the mean location of the two geoidal heights. The estimated standard deviation of g is simply
g d obtain geoidal height from SSH, both the time-dependent and quasi-time independent
SST values should be removed. In this paper, the time-dependent SST is reduced by filtering, and for the quasi time-independent SST we adopt the model of Levitus et al.
(1997), which is available at NOAA’s Ocean Climate Laboratory (see http://www.nodc.noaa.gov/OC5/dyn.html). The Levitus SST values from NOAA are given as monthly averages on a 1º×1º grid for 12 months. We averaged the monthly values to get the needed quasi time-independent SST. Figure 3.1 shows the averaged Levitus SST, which looks very similar to the 1982 version of Levitus SST (cf: Hwang, 1997, Fig. 3). The divided difference method (Gerald and Wheatley, 1994) is then used to interpolate the needed SST value from this SST grid at any altimeter data point.
Figure 3.1: Quasi time-independent sea surface topography from Levitus et al. (1997), contour interval is 10 cm.
3.2.2 Removing Outliers and Gridding DOV
Because of the use of 1D FFT algorithm in the computations of MSSH and gravity anomaly (see below), the north and east components of DOV need to be formed on regular grids with constant intervals in latitude and longitude. Before
forming the regular grids, erroneous along-track DOV must be removed. Again we adopt Pope’s tau-test method to remove possible outliers. Within a 4′×4′ cell, any along-track DOV, εi, can form an observation equation as (cf: Heiskanen and Moritz, 1985, p. 187)
ξ and are the north and east components. After least-squares estimating ξandη, all residuals in the 4′×4′ cell can be determined using Equation (3.4). A DOV is flagged as an outlier and is removed if its residual satisfies the condition:
) tau-value at degree of freedom of m; see Koch (1987, p.336) and Pope (1976) for the methods of computing σviand τc(m). Table 3.1 shows the ratios between removed and raw DOV in 12 selected areas. The removed DOV are largely from the non-repeat missions. In general, the removal ratios are relatively high in shallow waters and in higher latitudes, e.g., Ross Sea, and low in the open ocean, e.g., the Reykjanes Ridge and the East Pacific Rise. One problem with the above procedure of outlier removal is that, in areas with sparse data such as coastal regions and polar regions, there are not enough data points to produce large degree of freedom to make the result of the tau-test reliable, leading to undetected/improperly detected outliers.
After removing outliers the north and east components are then computed on a regular grid by the method of LSC (Moritz, 1980):
l
Table 3.1: Test areas and ratios of removed outliers
Area Geographic boundaries
where vector l contains along-tack DOV, vector s contains north (ξ) and east (η)components, Csl, Cll and Cn are the covariance matrices for s and l, l and l, and the noise of l, respectively. In Equation (3.6), Cn is a diagonal matrix that in theory contains the variances of along-track DOV. When a reference gravity model is used, the error of the gravity model must be taken into account in constructing the covariance functions, see Hwang and Parsons (1995) and Hwang et al. (1998) for the methods of constructing the covariance functions when gridding DOV by LSC.
3.3 Conversions from DOV to MSSH and gravity anomaly
We use the deflection-geoid and inverse Vening Meinesz formulae as the basic tool for computing MSSH and gravity anomaly. The deflection-geoid formula transforms DOV into geoidal height, which then yields MSSH by adding the quasi time-independent SST. Detailed derivations of these two formulae are given in Hwang (1998). These two formulae read
q
R: mean earth radius, 6371000 m is used
γ : normal gravity, based on GRS80 (Torge, 1989) H
C′and ′: Kernel functions :
, q
q η
ξ north and east components of DOV at q (dummy index)
αqp: azimuth from q to p σ : unit sphere
dσq: surface element = cosφqdφqdλ and φq,λqare latitude and longitude
The kernel functions C′andH′are functions of spherical distance only and are defined in Hwang (1998). The 1D FFT algorithm is used to rigorously implement
Equation (3.7). In the case of using a 360-degree reference field (see below), an optimal effective radius of integration in Equation (3.7) is about 110 km (about 1º at the equator). In the 1D FFT algorithm, all geoidal heights or gravity anomalies at a fixed latitude (or parallel) are computed simultaneously (Hwang, 1998), and this is why the 1D FFT algorithm is faster than the strait sum algorithm.
Because of the singularity of the kernel function C′andH′ at zero spherical distance, the inner most zone effects on geoidal height and gravity anomaly must be taken into account and are computed by
⎭⎬ s0 is the size of innermost zone, which can be estimated from the grid intervals as
π y s ∆x∆
0 = (3.9)
Formulae such as those in Equation (3.9) are based on spherical approximation.
Errors arising from spherical approximation are investigated in detail by Moritz (1980). When using the remove-restore procedure, error in using spherical approximations should be very small compared to data noises. Consider the formula of error-free LSC in the case of using ellipsoidal correction (Moritz, 1980, p. 328)
where s, l,C and sl Cll−1 are defined in Equation (3.6) and e2 is the squared eccentricity of a reference ellipsoid, which is about 0.006694 for the GRS80 ellipsoid.
The procedure in Equation (3.10) is first to remove the ellipsoidal effect of the data, l1
e2 , perform LSC computation and finally add back the ellipsoidal effect of the signal e2s1. In the case of using remove-restore procedure where a reference field is removed from the data and the residual signal is to be recovered, both l1ands1 will be very small compared to their full signals (see Moritz (1980), p. 327, where the low degree part will vanish due to the use of a reference field). For example, if the largest element (DOV) in l1 is 100 µrad, then the largest element in e2l1 will be 0.66 µrad. This value is far smaller than the noise of DOV from the multi-satellite altimetry. Furthermore, if the largest element in s1 is 100 mgal, then the largest ellipsoidal effect on gravity anomaly is 0.66 mgal, which is much smaller than the error of the recovered gravity anomaly. Of course, whether the reference field used in the remove-restore procedure will introduce additional error is another issue.
3.4 Computation and Analysis of Global MSS Model
After the tests performed in the previous sections and the selection of a set of optimal parameters, MSSH and gravity anomalies on a 2´×2´ grid were computed over the area: 80ºS-80ºN and 0º-360ºE. The computations were divided into 36 areas, each covering a 40º×40º area. A batch job was created for each of the 36 areas and this batch job creates maps of altimeter data distribution, predicted MSSH and gravity anomaly, and many statistics for detailed examinations. The final global MSSH and gravity anomaly grids are the combination of the results from the 36 areas. Figure 3.2 shows the flowchart of computation in a 40º×40º area. The most time-consuming part
in the procedure is the gridding of DOV. Only few minutes of CPU time is needed for the 1D FFT computation of geoid or gravity anomaly in a 40º×40º area on a Pentium III 600 machine. This computational procedure uses the EGM96 gravity model (Lemoine et al., 1998) to harmonic degree 360 as the reference field. In summary, reference DOV implied by EGM96 were removed from the raw DOV to yield residual DOV, which were used to compute residual geoidal heights and gravity anomalies using Equation (3.7). The final geoidal heights and gravity anomalies are obtained by adding back the EGM96-implied values.
As seen in Figure 3.2, the procedure for obtaining the global MSSH grid is more involved than that for the gravity anomaly grid. A preliminary MSSH grid was first obtained by adding the 1994 Levitus SST to the geoid grid. Because of the use of DOV as data type, it is possible that the long wavelength part of MSSH is lost in using the deflection-geoid formula. To mitigate such a loss, we first computed the differences between the along-track T/P, ERS-1 MSSH and the preliminary MSSH.
Each difference is associated with a weight, which is the inverse of noise variance.
Then smoothing and de-aliasing of the differences were made by computed the weighted median values within 15΄×15΄ cells. The weighted median values were then interpolated on a 15΄×15΄ grid using the minimum curvature method (Smith and Wessel, 1990). The final MSSH grid is obtained by summing the difference grid (which is now re-sampled into a 2΄×2΄ grid) and the preliminary grid. The resulting grids are now designated as the NCTU01 MSSH grid and the NCTU01 gravity grid.
Because of the final adjustment using T/P MSSH, The SSH from the global MSSH grid is the height above a geocentric ellipsoid with a semi-major axis=
6378136.3 m, and flattening=1/298.257222101. The geodetic coordinates of both
In addition, the normal gravity for the global gravity anomaly grid is GRS80 because we have removed the zonal spherical harmonic coefficients C20, C40, C60 and C80 of the GRS80 reference ellipsoid (cf: Torge, 1989) when computing the reference gravity anomalies from EGM96
A standard method for evaluating a MSSH grid is to compare modeled MSSH and MSSH-derived gradients with averaged SSH and gradients from repeat missions.
The first comparison is for the MSSH values from our model and from repeat missions. Table 3.2 shows the result of the comparison of MSSH values using the averaged along-track SSH from T/P and ERS-1. Also included in Table 3.2 are the comparisons for NASA/GSFC model (Wang, 2000) and the CLS model (Hernandez and Schaeffer, 2000). Compared to the NASA /GSFC and CLS models, the NCTU01 model agrees best with the T/P and the ERS-1 MSSH. The GSFC MSSH is slightly worse than the NCTU01 MSSH. On the continental shelves (depths below 200 m), all MSSH models contain large errors, which in the case of CLS has exceeded 20 cm.
Even in the median depths, the accuracy of MSSH is not very promising and is generally worse than 10 cm. All MSSH models have the best accuracy of few cm in the deep oceans. The differences in the case of ERS-1 are smaller than in the case of T/P. This is to be explained by the fact that ERS-1 has a higher data density than T/P, so the former will dominate the resulting MSSH model. As such, the MSSH model will have a better match with the ERS-1 SSH than T/P SSH. Furthermore, Figure 3.3 shows the differences between the NCTU01 and T/P MSSH. Again the differences in Figure 3.3 are large over shallow waters and small in the deep oceans. It is noted that the pattern of differences is very similar to the pattern of ocean variability in Figure 2.3.
remove DOV of
Seasat SSH Geosat SSH ERS-1 SSH ERS-2 SSH T/P SSH
north and east DOV
MSSH grid gravity anomaly grid
Figure 3.2: Flowchart for computing global MSSH and gravity anomaly grids
This shows that the MSSH is less reliably determined in areas of high ocean variability than in other areas. Such a result agrees with the expected outcome of least-squares collocation (see Equation (3.6)): data with larger noises (variability) yield less reliable results.
Table 3.2: RMS differences (in cm) between global sea surface models and T/P and ERS-1 MSSH
T/P ERS-1
0-0.2a 0.2-1 1-10 0-10 0-0.2 0.2-1 1-10 0-10
NCTU01 16.5 10.6 1.9 5.0 6.8 4.2 2.7 3.1
CLS 24.7 26.1 3.9 9.0 12.2 8.8 4.7 5.4
GSFC00 19.4 13.5 2.0 6.0 10.7 7.6 3.6 4.3
No. points 34824 23015 485762 543601 63828 48296 1517901 1630025
a Depths from 0 to 0.2 km
Figure 3.3: Difference between NCTU01 and T/P mean sea surface heights
Table 3.3 shows the RMS differences in gradients between the NCTU01, GSFC and CLS modeled gradients and the averaged along-track gradients from T/P and ERS-1 missions. The NCTU01 and GSFC models have almost the same level of accuracy, while the CLS model has a slightly poorer accuracy than the two. For all models the distribution of error is the same as that of MSSH error, namely, bad accuracy over shallow waters and good accuracy in the deep oceans. Again, for the reason of data dominance, the error in the ERS-1 case is smaller than the error in the T/P case. In the T/P case the RMS difference of 7-9 µrad has far exceeded than the standard deviation of T/P gradient, which is 1.9 µrad based on Equation (3.3) (with an averaged point spacing of 6.6 km, see Table 2.2). On the other hand, in the ERS-1 case the RMS difference of 3.2 µrad agrees well with the expected standard deviation of 4.7µrad of ERS-1 gradient.
Table 3.3: RMS differences in along-track sea surface height gradients (in µrad) derived from global sea surface models and from T/P and ERS-1 MSSH
T/P ERS-1
0-0.2a 0.2-1 1-10 0-10 0-0.2 0.2-1 1-10 0-10
NCTU01 23.3 16.7 2.4 7.1 7.1 5.0 2.8 3.2
CLS 26.4 19.9 2.9 9.0 8.7 5.3 2.8 3.3
GSFC00 22.9 16.9 2.3 7.0 7.3 4.3 2.8 3.2
No. points 33687 22814 485561 542062 61374 47071 1510808 1619253
a Depths from 0 to 0.2 km
3.5 Computation and Analysis of Global Gravity Anomaly Model
To evaluate the accuracy of the NCTU01 gravity grid, we made comparison between the predicted and shipborne gravity anomalies in 12 areas in the world oceans (see Table 3.1 for the boundaries). These shipborne gravity data are from the National Geophysical Data Center. The corresponding ship tracks are shown in Figure 3.4. The 12 areas are so selected that deep and shallow waters, and low and high latitudes are covered in the comparison. The comparing procedure and the adjustment of shipborne gravity data were detailed in Hwang et al. (1998). Briefly, long wavelength biases in the shipborne gravity anomalies were first adjusted using altimetry-derived gravity anomalies. A pointwise comparison was then made for the adjusted shipborne gravity anomalies. Table 3.4 shows the result of the comparison.
In Table 3.4, the gravity grid from Hwang et al. (1998) is also compared. It has been shown that the gravity grid of Hwang et al. (1998) has a better accuracy than that of Sandwell and Smith’s (1997). As shown in Table 3.4, the accuracy of the current gravity grid has indeed improved over the gravity grid of Hwang et al. (1998). The improvement is most dramatic in the deep oceans such as the East Pacific Rise. In the Caribbean Sea, the improvement is also significant. Other areas with significant improvements are the Ross Sea and Kerguelen Plateau, which are situated in the Antarctica. In general, large difference between predicted and shipborne gravity anomalies may be due to either or both of the following reasons: (1) bad altimeter
In Table 3.4, the gravity grid from Hwang et al. (1998) is also compared. It has been shown that the gravity grid of Hwang et al. (1998) has a better accuracy than that of Sandwell and Smith’s (1997). As shown in Table 3.4, the accuracy of the current gravity grid has indeed improved over the gravity grid of Hwang et al. (1998). The improvement is most dramatic in the deep oceans such as the East Pacific Rise. In the Caribbean Sea, the improvement is also significant. Other areas with significant improvements are the Ross Sea and Kerguelen Plateau, which are situated in the Antarctica. In general, large difference between predicted and shipborne gravity anomalies may be due to either or both of the following reasons: (1) bad altimeter