Zonal and Meridional Ocean Currents at TOPEX / Poseidon and JASON-1 Crossovers around Taiwan : Error Analysis and Limitation

A crossover method for determining zonal and meridional ocean current components is examined using data at three crossovers of TOPEX/Poseidon and JASON-1 ground tracks over 2002 2006. To implement this method, a geoid model around Taiwan is constructed using surface and airborne gravity data. The modeled and observed geoidal heights at coastal benchmarks are consistent to 5 cm RMS with the means removed. The error and limitation of this method are discussed, concluding that, in order to obtain current velocities at a 10 cm s-1 accuracy and a 6-km resolution, the dynamic ocean topography (DOT) at a mmlevel accuracy is needed, which is not possible to achieve today. By filtering DOT to a spatial scale of 100 km or coarser, a 10 cm s-1 accuracy of velocity may be obtained. One crossover (A) is situated south of Taiwan and near the Kuroshio, the second (B) is at the axis of the Kuroshio and the third is located in the northern Taiwan Strait. These three crossovers feature different ocean current patterns. At a spatial scale of 120 km, the agreement among the altimeter, the Princeton Ocean Model (POM), and the drifter-derived velocities is the best at B, followed by that at A, and then C. In fact, at C the altimeter-derived velocities contradict the POM-derived values, and the tide model error is to be blamed. Further improvement on geoid modeling is suggested.


INTrOduCTION
Tools for observing ocean currents can be surface and space-based.One of the space-based tools is satellite altimetry, which has not only delivered a revolutionary result in oceanography, but also new findings in geophysics and geodesy.A comprehensive document of techniques and applications of satellite altimetry can be found in Fu and Cazenave (2001).Oceanographic studies around Taiwan based on altimetry are abundant in the literature.The most frequently investigated subjects, among others, are the Kuroshio and eddies east of Taiwan, tides, circulations and eddies in the South China Sea, and tides in the Taiwan Strait.For example, 18 papers are found in the Web of Science Database (http://portal.isiknowledge.com) using the keywords: Taiwan, Kuroshio, and altimetry.Due to shallow waters and bad altimetry data quality, ocean currents in the Taiwan Strait solely derived from altimetry are, however, not seen in the literature.
Varying at different spatial and temporal scales, ocean currents around Taiwan are rather diversified (Liang et al. 2003).Originating in the tropical Pacific, the Kuroshio lies east of Taiwan and is one of the major western boundary currents of the world's oceans.The South China Sea is situated south of Taiwan and is a semi-closed sea, and here the circulations are largely affected by monsoonal winds and are therefore mostly seasonal.Both the Kuroshio area and the South China Sea are abundant with warm and cold-cored eddies (Hwang et al. 2004).The Taiwan Strait lies west of Taiwan and is a newly explored area in terms of oceanography and here numerous findings in the 2000s have filled information gaps; see, e.g., Chen and Sheu (2006).North of Taiwan lies the East China Sea, where the Kuroshio "in-trudes" into the continental shelf in the winter time (Tang et al. 2000).These four oceanic zones are now under intensive studies in the world oceanographic community.
This paper does not attempt to study all these ocean currents by satellite altimetry.Rather, we will investigate an altimetry-based method (called the crossover method) that can resolve the along-track dynamic ocean topography (DOT) into the zonal and meridional ocean current components at the crossover of two satellite ground tracks.The potential of this method to transform altimeter observations into ocean current observations at a fixed point and at a regular time interval is examined.We choose to use TOPEX/ Poseidon (T/P) and JASON-1 altimeter data.These two missions are specifically designed for oceanographic applications.The 10-day repeat period is ideal for seeing shortterm variations of oceanic signals.Our method requires an accurate geoid model and accurate sea surface heights (SSHs) from altimetry.In the following development, we will show the gravity data and the method of geoid modeling around Taiwan.The error and limitation associated with the crossover method of ocean current determination will be discussed.

ZONAL ANd MErIdIONAL VELOCITy COM-PONENTS AT CrOSSOVEr POINTS
The ocean current investigated in this paper is restricted to the geostrophic current, which is a result of the balance between the pressure gradient and the Coriolis force (Cushman-Roisin 1994).In this case, the zonal and meridional velocity components are determined by: where g is gravity at the crossover point, sin with ω e being the mean rotational velocity of the earth and φ latitude; x and y are the zonal and meridional coordinate variables; and ξ is DOT.A point DOT can be determined by: where h is sea surface height from satellite altimetry, and N is geoidal height.At the crossover of two satellite ground tracks, the gradients of DOT along the ascending and descending tracks, ε a and ε d , can be expressed by the east and north gradient components as (Fig. 1; Inverting ( 4) and (5) yields: which can be used in Eqs. ( 2) and ( 3) to compute the velocity components.This is called the crossover method.The formula of Parke et al. (1987) that decomposes the alongtrack velocities into orthogonal components is similar to Eqs. ( 6) and ( 7).The Parke et al. (1987) formula is based on an approximate relationship between the azimuths of the ascending and descending satellite ground tracks.Morrow et al. (1994) and Strub et al. (1997) use the same formula as that given in Parke et al. (1987) for their studies.
A crossover may come from a single and a dual satellite mission.Use of altimetry data from a single satellite mission yields a velocity at a crossover of two satellite ground tracks at a time interval equal to the satellite repeat period.Therefore, the T/P, Geosat-Follow-On (GFO) and ENVI-SAT missions will yield temporal resolutions of ocean current at 10, 17, and 35 days, respectively.It is also possible to use crossovers from a dual satellite mission, but the temporal resolution will be degraded to that of the satellite mission with the longest repeat period.For example, use of T/P-ENVISAT crossovers will yield a temporal resolution of 35 days.Differences in the system accuracies, inclinations and spatial resolutions of the involved satellite missions may produce inconsistent accuracies in the ocean current components.Since T/P and JASON-1 have the same repeat period and inclination, the temporal resolution of ocean currents derived from these two missions remains 10 days.

ErrOr ANALySIS ANd ACCurACy ChAL-LENgE
It turns out the accuracy requirement for DOT in the crossover method is very demanding, as given in a theory below.In order to see the relationship between the error of DOT and the error of geostrophic velocity, we approximate an along-track gradient of DOT by the finite difference: .
where Δh, ΔN, and Δs are differenced SSH, geoidal height and along-track distance between two consecutive points, respectively.A differenced SSH is less affected by long wavelength errors than the original SSHs (Hwang and Parsons 1996).Examples of long wavelength errors are the one cycle per revolution (CPR) orbit error and the tide model error in the open ocean.If there is no correlation between Δh and ΔN, the standard error of gradient is: We assume that the standard errors of the gradients of DOT for the ascending and descending tracks are the same.By the approximation (Sandwell 1992): the standard errors of the zonal and meridional velocity components can be determined as: where α = α a is the azimuth of the ascending track.Thus, the ratio between the two standard errors is: Using the theory of azimuth given in Hwang and Parsons (1996), one can derive an approximate expression of azimuth as: where , u X o o are the angular velocities of the argument of latitude and the ascending node of the satellite in the inertial frame, I is the inclination of satellite orbit (T/P and JASON-1 are about 66°).If the orbit undergoes a linear perturbation caused by the degree 2 zonal harmonic, J 2 , then , u X o o and I will be nearly constants (Kaula 1966).Therefore, the errors of velocity components [Eqs.( 11), ( 12)] and the ratio r are functions of latitude.The phenomenon that the velocity error depends on latitude was empirically stated in Schlax and Chelton (2003), and it is now explained by the above formulae.Figure 2 shows the relationships between r and latitude for the T/P (also JASON-1), ENVISAT and GFO missions.For these missions, the error of the zonal component is always smaller than that of the meridional component.
While r is independent of spatial scale, the errors of the velocity components are not.In Eqs. ( 11) and ( 12), Δs works as a filtering parameter.At a mean latitude of 23.5° (this is the case for Taiwan), we have 22 a = c for T/P.Assuming 2 (Fu et al. 1994) and  5 N v = D cm (see Section 5 below), the standard errors of velocity components follow the empirical relationships: . s . s where the standard errors are in ms -1 , and Δs is in km.For the one HZ T/P altimeter data, we have Δs = 6.2 km, leading to the result: . 1 18 These errors exceed the magnitudes of most ocean current velocities in the world' oceans.On the other hand, to achieve a 10 cm s -1 accuracy in velocity, we will need Δs = 73 km for the zonal component and Δs = 207 km for the meridional component.However, if we can achieve a mm-level accuracy for DOT from satellite altimetry and a geoid model, then we can achieve a 10 cm s -1 accuracy in velocity at spatial scales of 7 to 20 km.Given the current altimetry technology and geoid modeling technique, it is not likely to obtain a mm-level accuracy for DOT.Therefore, only the mesoscale geostrophic velocity components (roughly at a spatial scale of 100 km or larger) at a 10 cm s -1 accuracy can be achieved using the crossover method.
The above error analysis can be extended to the case of a two-dimensional geostrophic velocity field determined using Eqs.( 1) and ( 2).Assume that the north and east components of the DOT gradient have the same accuracy.Applying error propagation to Eqs. ( 1) and ( 2) yields the following standard error of velocity: Again, for 2 cm, and Δs = 6.2 km, we get . 1 58 m s -1 at a mean latitude of 23.5°.In order to achieve, for example, 1 u v v v = = and 10 cm s -1 , the spatial scales will have to be Δs = 1000 and 100 km, respectively (latitude = 23.5°).The case Δs = 100 km corresponds to a degree 360 field in the spherical harmonic expansion of the geopotential, or the geoid.Recent geopotential models from such state-of-the-art mission as GRACE (Tapley et al. 2005) still cannot deliver a geoid model accurate to 5 cm at this wavelength (100 km).

TOPEX/POSEIdON ANd JASON-1 ALTIMETEr dATA ANd CrOSSOVErS
The altimeter data for estimating the ocean current velocities are from the T/P and JASON-1 missions, and are provided by AVISO (1996); see also http://www.aviso.oceanobs.com.The first JASON-1 repeat cycle corresponds to T/P cycle 344.Starting with Cycle 365, T/P was maneuvered into a new orbit whose ground track is between the original and current JASONS-1 ground tracks.The new T/P repeat mission began with Cycle 369 (September 20 2002).In this paper we use data from T/P cycles 367 to 480 and JASON-1 cycles 19 to 138.Our altimeter data processing uses standard geophysical corrections for the raw T/P and JASON-1 SSHs, namely, tide (CSR4.0tide model), wet troposphere (from the onboard radiometer result), dry tropo-sphere (ECMWF), ionosphere, inverse barometer, sea state bias, pole tide, and drift corrections.
Figure 3 shows the distribution of T/P and JASON-1 crossovers and the bathymetry around Taiwan.Based on the availability and the accuracy of our geoid model, we choose crossovers A, B, and C to compute velocities.Table 1 show the facts about A, B, and C. A and C are crossovers of T/P and JASON-1, and B is a crossover of T/P only.These three crossovers are more than 30 km to the coasts, so the altimeter ranging will not be degraded by waveform corruption; see also Deng (2004).C is over shallow waters, so here the SSHs may contain a larger tide model error than those at A and B. According to Liang et al. (2003), A and B are near the Kuroshio, where the ocean current is largely geostrophically balanced and the maximum speed may exceed 100 cm s -1 .A is located in a zone where the Kuroshio meanders.C is in the northern end of the Taiwan Strait Current (Liang et al. 2003;Wu et al. 2007).
Other crossovers in Fig. 2, not explored in this paper, may be useful for ocean current research around Taiwan.For example, the ocean current at D can be used to investigate the seasonal variation of the Kuroshio northeast of Taiwan (Zhang et al. 2001;Hwang and Kao 2002).The ocean current at E is useful to see the exchange of water between the Kuroshio and the South China Sea through the Luzon Strait Fig. 3. Distribution of T/P and JASON-1 crossovers around Taiwan.(Caruso et al. 2006;Wu et al. 2007).F is an ideal location to study the circulation of the northern South China Sea, and the ocean current here may be used to validate or disprove the entrance of Kuroshio water into the Taiwan Strait.

grAVITy dATA
As stated before, our method of ocean current determination requires an accurate geoid model.To this end, we collected all possible gravity data around Taiwan.The gravity data contain surface, airborne and altimeter-derived data.The surface gravity data include land and shipborne data (Fig. 4a).Land gravity data were collected over 1980 -2004 by several agencies, including Academia Sinica, the Base Survey Battalion and Ministry of Interior (MOI) of Taiwan.The standard errors of land gravity range from 0.04 to 1 mgal (Yen et al. 1995;Hwang et al. 2002).The shipborne gravity data were collected by international research vessels from 1960s to 1990s, and these gravity data are available at the National Geophysical Data Center (NGDC; http:// www.ngdc.noaa.gov).The shipborne gravity data we used are from NGDC and have been edited by Hsu et al. (1998).The qualities of these data vary from one cruise to another.In general, measurements collected in 1990s have a better accuracy than those in the 1960s to 1980s.For example, gravity data collected by Taiwan's Ocean Researcher 1 vessel in 1996 using an Atalante KSS30 gravimeter have a mean and standard deviation of 0.3 and 2.6 mgal in the crossover differences of gravity values.According to Hsu et al. (1998), the mean and standard deviations of crossover differences of all the shipborne gravity data are 6.3 and 11.2 mgal, respectively.
The airborne gravity data were from a survey carried out in 2004 and 2005 (Hwang et al. 2007).The numbers of survey lines, flight altitude and speed of aircraft are 100, 5156 m and 300 km hr -1 , respectively.The cross-line intervals are 4.5 and 20 km for the north-and east-going lines, respectively.A LaCoste and Romberg System II Air-Sea gravimeter was used to collect the data.The average RMS crossover difference of gravity at intersecting lines is about 2.9-mgal, equivalent to a 2-mgal point standard error.Due to filtering, the spatial resolution of the airborne gravity data is about 6 km.Finally, the altimeter-derived gravity anomalies (Fig. 4c) were from Andersen et al. (1999) and were derived from the Geosat/GM and ERS-1/GM.A 5-mgal standard error was assigned to all altimeter-derived gravity anomalies.

CONSTruCTINg ANd VALIdATINg A gEOId MOdEL ArOuNd TAIwAN
The method for geoid modeling is based on the remove-compute-restore (RCR) procedure.The computation of a geoidal height, N, is divided into three parts: where N ref is geoidal height from a geopotential model, N res is the residual geoidal height and N rtm is geoidal height due to the residual terrain model (RTM) (Forsberg 1984).Accordingly, the surface (including altimeter-derived) and airborne gravity anomalies Δg surf and Δg air can be also spilt into three portions:   For the reference geopotential model, we adopt a combined GGM02C (Tapley et al. 2005) and EGM96 (Lemoine et al. 1998) model to degree 360.The elements of the covariance matrices in Eqs. ( 21) and ( 22) are determined from covariance functions at given spherical distances of two data points.The covariance functions are expanded into series of Legendre polynomials.For the covariance function of the disturbing potential, the coefficients of expansion are the error degree variances of the combined GGM02C-EGM06 model (for degree 2 to 360) and the model degree variances of Tscherning and Rapp (1974, Model 4, for degree 361 to infinity).Other covariance functions are derived by covariance propagation, see also Hwang and Parsons (1995).For the RTM effects, we use a 3" × 3" digital elevation model of Taiwan.Figure 5 shows the geoid model and the error distribution.The geoid variation and the formal errors are correlated with the roughness of terrain and bathymetry.The formal errors at A, B, and C are 34, 22, and 5 cm, respectively.This difference in formal error is partly due to the fact that A and B are over a rough bathymetry and C is over a smooth bathymetry, and partly due to different densities of gravity data.Note that the airborne gravity data contribute mostly to the geoidal signal at C. It is rather difficult to carry out an external accuracy assessment of a geoid model at sea.As a compromise, we did an assessment on the benchmarks along two leveling routes in the north coastal and east coastal zones close to A, B, and C.Here the modeled and the observed geoidal heights are compared, the later being defined as the difference between the GPS-derived ellipsoidal height and the orthometric height from precision leveling.Table 2 shows the statistics of the comparison.The standard deviation of the differences in Table 2 is a descriptor of the geoid accuracy.Again, due to the smooth bathymetry and smooth gravity variation, the geoidal heights near C (north route) are better determined than those near A and B (east route).According to Table 2, we could expect a geoidal accuracy of about 5 cm near A, B, and C.

rESuLT Of VELOCITy dETErMINATION ANd VALIdATION
Using the DOT from the T/P and JASON-1 altimeter data and the geoid model (Section 6), we have computed zonal and meridional velocity components at crossovers A, B, and C at a 10-day interval.According to the error analysis in Section 3, it is not possible to obtain a reliable velocity using un-filtered DOTs.Therefore, we used the Gaussian filter of GMT (Wessel and Smith 1999; http://gmt.soest.hawaii.edu/) to do the filtering.A filtered DOT is the result of the convolution of the Gaussian function / ) s d ( exp 2 2 -with raw DOTs within a window of filter width that is equal to Due to missing satellite passes and/or missing height observables near the crossovers, these time series all contain gaps.The altimeter-derived velocities were then compared with in-situ velocities computed from drifter data.The drifter data are from the World Ocean Circulation Experiment (WOCE) database (http://woce.nodc.noaa.gov/wdiu/).The drifter data are rather sparse around A and B, and there is no data around C (Fig. 9).A drifter-derived velocity was computed as the ratio of traveling distance and the difference of times between two successive records.Because a drifter is not likely to exactly pass through any of the crossovers, we use drifter data within a 0.1° box centered at the crossover.Note that the drifter velocities may be sub-ject to systematic errors and random noises caused by tidal currents and positioning and timing errors.In general, the difference between the altimeter-derived and the drifter-derived velocities decreases as the filter width increases.The altimeter-derived velocities were also compared with the model-derived velocities of Wu and Hsin (2005).The model of Wu and Hsin (2005), driven by winds and sea surface temperature, is based on the Princeton Ocean Model (POM; Blumberg and Mellor 1987) and the boundary conditions of East Asia Marginal Seas (EAMS).Below is a summary and discussion based on the results in Figs. 6 to 8.

Crossover A
At a filter width < 120 km, the altimeter-derived velocities are larger than the model and drifter-derived velocities, and the directions of flows from three results are reasonably consistent.At filter width = 120 km, the magnitudes of velocities from all results differ by about 10 cm s -1 .At any filter width, the directions of flows from all results are mostly northward or northeastward.In rare cases the altimeter and model-derived currents turn southwards.The drifterderived currents are nearly northward in most cases.Since the Kuroshio meanders around A, the ocean currents here fluctuate rapidly (C.W. Wu, private communication, 2006).This fluctuation is also seen in Fig. 6.   2).

Crossover b
Here in most cases the altimeter-derived velocities exceed 200 cm s -1 and are larger than the drifter and model-derived velocities at any filter width.In some cases, the drifterderived velocities are nearly 200 cm s -1 , otherwise fall be-tween the altimeter and model-derived velocities.Compared to A, the velocities at B from the altimetry and the POM model appear to be more time independent, and the directions of flows are uniform throughout the entire year.Note that the altimeter-derived velocities (filter width = 120 km) at A and B agree well with the velocities given in Liang et al. (2003).error over the Taiwan Strait should be one reason for the discrepancy between the altimetry and the model results.For example, the crossover differences of along-track SSHs at C are found to be as large as 60 cm, which should be mainly caused by tide model error.For comparison, the crossover differences at A and B are less than 10 cm.Also, C is about

Crossover C
At the filter widths of 40, 80, and 120 km, the altimeter-derived velocities at C do not match the model-derived velocities in both magnitudes and directions.The altimeterderived currents are mostly southeastward, while the modeled currents are mostly northeastward.A large tide model ocean current velocities here is to assimilate altimeter data into such ocean dynamic models as POM; see also the discussion in Le Provost and Bremond (2003).
Reasonably good ocean currents at a spatial scale larger than 100 km were obtained from T/P and JASON-1 altimetry at crossovers A and B (the Kuroshio area), but the 31 km to the coast and here the depth is 79 m, so the condition of geostrophic balance, which is valid for the open ocean, may not hold here.That is, the method in Section 2 may not be valid at C. Furthermore, we cannot blame waveform corruption near C for the inconsistency between the altimetry and the model results, because C is more than 10 km away from the coast (Deng 2004).

Interannual variation
The time series in Figs. 6 to 8 have been spatially filtered.In order to see the temporally filtered velocities and in particular the interannual variation, we applied the same Gaussian filter to these time series using a filter width of 12 months for the 120-km filtered time series.The results are given in Fig. 10.The time-filtered, altimeter-derived velocities seem to be rather smooth at all crossovers.At A and C, the pattern of velocity did not repeat exactly from one year to another.At A, the ocean currents in 2005 were the weakest.Note that a La Niña occurred over 2005 -2006.An interannual variation of velocity is also evident at C and again the currents in 2005 here were also the weakest.At B, the ocean currents are quite stable interannually.A possible explanation for this stability is that, B is almost at the axis of the Kuroshio, which maintains a steady flow from one year to another.

dISCuSSION ANd CONCLuSION
This paper presents a crossover method for computing zonal and meridional ocean current components.An error analysis leads to the conclusion that, in order to obtain ocean current velocities at a 10 cm s -1 accuracy and a 6-km resolution, we need a DOT at a mm-level accuracy.Given the current altimetry technology and the accuracy of geoid modeling, to achieve a 10 cm s -1 accuracy for velocity the DOT must be filtered to a 100-km spatial resolution or coarser.This implies that only mesoscale ocean currents can be obtained by the crossover method (Section 3).So the limitation of the crossover method is on the spatial scale.Note the accuracy of such altimeter-velocity is also latitude dependent.
The geoid model obtained in this paper can also be used to compute a two-dimensional geostrophic velocity field around Taiwan.That is, Eq. ( 3) can be used to compute a grid of DOT using altimeter data from multi-missions such as T/P, JASON-1, ERS, and ENVISAT.Even data from the laser altimeter mission ICESat (Schutz et al. 2005;Urban et al. 2008) can be used to enhance data coverage, especially around coasts.Again, given the current DOT accuracy, the spatial resolution of a 2-D velocity field will be at the mesoscale level (coarser than 100 km).Since the condition of geostrophic balance does not strictly hold over shallow waters (e.g., the Taiwan Strait), a better approach of obtaining The legend is the same as that in Fig. 6. result at C (the Taiwan Strait) contradicts the model output of POM (Section 7).It is possible to improve our results at A, B, and C by taking into account the following points: 1) The marine gravity data must be improved in both coverage and accuracy.The systematic and random errors associated with the marine gravity data used in this paper must be further calibrated; 2) The current technique of geoid modeling is optimal for land, but might not be so for the sea.An improved marine geoid around Taiwan is needed; 3) An improved tide model over the Taiwan Strait is needed.This can be achieved by assimilating tidal constants at tide gauges and altimeter SSHs into a hydrodynamic model with a fine-resolution bathymetry model.

Fig. 1 .
Fig. 1.Geometry showing the crossover of two satellite ground tracks and velocity components.

Fig. 2 .
Fig.2.Ratio between errors of zonal and meridional velocities as a function of latitude for the T/P, ENVISAT and GFO satellite missions.
at sea level and the flight altitude.The residual geoid, N res , and the error variance were computed by: is the variance of geoid; C ng surf , C ng air , C g surf , C g g surf air , C g air , and C g g surf air are covariance matrices for geoid-surface gravity anomaly, geoid-airborne gravity anomaly, surface gravity anomaly -surface gravity anomaly, airborne gravity anomaly -surface gravity anomaly, airborne gravity anomaly -airborne gravity anomaly, surface gravity anomaly -airborne gravity anomaly; and D g surf D and D g air D are the error variances (see Section 5) of the surface and airborne gravity anomalies.
Fig. 10.Interannual variation of velocity at A (top), B and C (bottom).The legend is the same as that in Fig.6.

Table 1 .
Facts about crossovers A, B, and C.

Table 2 .
Statistics of differences (in meter) between observed and modeled geoidal heights at two leveling routes near coasts.