Calibrating error variance and scaling global covariance function of geoid gradients for optimal determinations of gravity anomaly and gravity gradient from altimetry

With the advancements of satellite altimeter technology and data volume, the accuracy and spatial resolution of altimeter-derived marine gravity fields have been steadily improved. Ideally, a best marine gravity field from multi-altimeter missions is one that combines all data by optimally calibrating their relative weights. To this end, we use the minimum norm quadratic unbiased estimator theory to calibrate error variances of geoid gradients (GGs) from the Cryosat-2 mission and the Jason-1/2 geodetic missions and then use a scaling factor to modify the global covariance functions into local covariance functions. The calibrated error variances and the scaling factors are used in the least-squares collocation to estimate the gridded north and east components of GGs, which are used to compute gravity anomaly and vertical gravity gradient (VGG) by the methods of inverse Vening-Meinesz (IVM) and numerical differentiations. The assessment of the altimeter-derived gravity anomalies in the northern part of the South China Sea using the shipborne gravity data shows an average gravity accuracy improvement of 9.5% by calibrated and scaled covariances of GGs compared to the initial variances. The method for VGG computations is confirmed by examining the extinct ridges in the Gulf of Mexico.


Introduction
Satellite altimetry has many applications in geodesy and earth sciences. As the altimeter technology advances, the accuracies and spatial resolutions of sea surface height (SSH) measurements improve over time. SSHs from satellite altimeters have been used to determine marine gravity anomalies since the era of Seasat in the 1970s (Haxby 1987;Hwang and Parsons 1995;Sandwell et al. 2014;Yu et al. 2021). Several global grids of marine gravity are now freely available, e.g., DTU17 (Andersen and Knudsen 2020) and Sandwell V31.1 (Sandwell and Smith 1997), vastly advancing marine geophysics and coastal gravity fields for geoid modeling.
In the 1990s, satellite altimetry underwent a substantial advance in data density and accuracy, which motivated many B Cheinway Hwang cheinway@nycu.edu.tw; cheinway@gmail.com 1 Department of Civil Engineering, National Yang Ming Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan, ROC researchers to develop improved methods for marine gravity derivations. The data and method improvements have led to a number of global marine gravity grids (Andersen andKnudsen 1996, 1998;Hwang et al. 1998Hwang et al. , 2002Hwang andParsons 1995, 1996;Sandwell and Smith 1997;Schwarz et al. 1990;Wang 1999Wang , 2001Andersen 2013). When combining altimeter data from multiple missions for gravity determination, proper weights for the individual mission data are important for achieving an optimal result. For example, Sandwell and Smith (1997) used a priori error estimates of the along-track slopes for weighting the data from Geosat and ERS-1 geodetic missions (GMs). Hwang et al. (1998Hwang et al. ( , 2002 used empirical variances for SSHs from the missions of Seasat, Geosat exact repeat mission (Geosat/ERM), Geosat/GM, ERS-1/ERM, ERS-1/GM, ERS-2, and TOPEX/POSEIDON data. Andersen and Knudsen (1998) created global, gridded SSHs of Geosat and ERS-1 GM data using the fast Fourier transform (FFT) without specifically stating the standard errors of the SSHs. The gravity recovery approach of Wang (2001) was similar to that of Andersen and Knudsen (1998) but was based on data from Geosat/ERM, Geosat/GM, ERS-1/ERM, ERS-1/GM, ERS-2, and TOPEX/POSEIDON without specific error estimates for the data.
Since 2010, new GM data from Cryosat-2, Jason-1/2, and SARAL/AltiKa have been available and were incorporated with the earlier altimeter data from Geosat and ERS-1 to recover global and local marine gravity grids with higher accuracies and resolutions than for previous grids. For example, Zhang et al. (2017) determined gravity anomalies over Southeastern China Seas from some of these recent altimeter datasets. Sandwell et al. (2019) derived global gravity fields from all existing altimeter GM data. Zhu et al. (2020) determined marine gravity over the South China Sea (SCS) from multi-satellite GM/ERM altimeter data and estimated the weights of geoid gradients (GGs) from the Ka-band SARAL/AltiKa altimeter and other Ku-band missions. In these sample studies, standard errors of the SSHs from an altimeter mission are not specifically given and are not calibrated with each other for best relative weights. The use of such uncalibrated weights may not lead to an optimal marine gravity field. Thus, optimizing data weights for altimeter observations from multiple missions can be an improvement for a marine gravity derivation method, especially when a statistical approach such as least-squares collocation (LSC) is used to construct grids of north and east gradient components, which can be used to determine marine gravity anomalies and vertical gravity gradients (VGGs).
The VGG is the vertical derivative of gravity anomaly, and it can amplify short-wavelength information in the gravity anomaly (Wessel and Lyons 1997). Both marine gravity anomalies and VGGs can be used to estimate oceanic depths (Smith and Sandwell 1994;Wang 2000;Hsiao et al. 2016), but the VGG is more sensitive than gravity anomaly to shallow anomalous masses and geologic structures. As such, the VGG is better suited to mapping relative fine-scale seafloor tectonic structures (Marson and Klingele 1993;Li 2001;Harper et al. 2021). The VGG has been used to detect the tectonic signatures associated with continental margins in the SCS  and the Gulf of Mexico (Sandwell et al. 2014). There are several methods for computing VGGs. For example, Rummel and Haagmans (1990) determined VGGs at crossover points from the second derivatives of altimetric geoid heights by Laplace's equation. Sandwell and Smith (1997) computed VGGs from the horizontal derivatives of the east and north vertical deflections by numerical differentiations. Wang (2001) calculated VGGs from gravity anomalies using the integral formula given in Heiskanen and Moritz (1985). Bouman et al. (2011) discussed possible systematic effects in the computation of VGG from altimeter data. Again, like previous studies for gravity derivation, these papers for VGG derivations used uncalibrated altimeter data noises. Bouman (2012) presented rigorous relationships between geoidal undulation, deflection of the vertical (DOV) and vertical gravity gradients.
In view of the potential benefit of improved altimeter data weights for gravity and VGG recoveries, the first objective of this paper is to use the minimum norm quadratic unbiased estimator (MINQUE) (Rao 1971;Koch 2013) to calibrate multi-satellite error variances. With ever-increasing altimeter data densities and redundancies for parameter estimation (see Sect. 3), the use of a method like MINQUE for altimeter noise calibration becomes possible. The second objective is to use a simplified method for scaling the global covariance functions of GG for a best LSC estimation of north and east gradient components, which are then used to determine marine gravity anomalies and VGGs in the Gulf of Mexico and the northern part of the SCS. The resulting gravity anomalies will be assessed by shipborne measurements, and the VGGs will be validated by the result of Sandwell V31.1 VGGs. When using the gridded north and east gradient components to compute VGGs, we also consider small terms that had been neglected by other authors.

Satellite altimeter data
To experiment with the methods of calibration, we use altimeter measurements from the Cryosat-2 and the Jason-1/2 GM. Cryosat-2 data are from the Baseline-D Geophysical Data Record (GDR) of the European Space Agency (ESA). Jason-1 data are from the GDR version "e" native products of the Centre National d'Etudes Spatiales (CNES). Jason-2 data are from the GDR version "d" products of CNES. Table 1 lists the information about these altimeter data and missions.

Resampling the 20-HZ SSHs to 2-HZ SSHs
To calibrate the relative weights of GGs from the three altimeter missions, we need initial standard errors of GGs, which are determined from the standard errors of SSHs. First, we resampled the 20-HZ SSH measurements from the three altimeter missions into 2-HZ SSH measurements and estimated their standard errors. The standard errors of SSHs were then used to estimate the standard errors of GGs by error propagation. The Cryosat-2 GDRs directly provide the 20-HZ SSHs with instrument, range, and geophysical corrections. The 20-HZ SSHs from Jason-1 and Jason-2 were obtained by subtracting the range measurements from the altitude measurements, and then applying the corrections recommended in the data manual, including range and geophysical corrections. The range corrections contain the wet/dry troposphere correction, ionosphere correction, and sea state bias correction. The geophysical corrections include tide effects (the sum of the solid earth tide, ocean tide, and pole tide), inverted barometer height correction, and high-frequency fluctuations of the sea surface topography.
To obtain smooth SSHs while retaining shortwavenumber information, we fitted 3-s SSHs (60 samples) by a quadratic polynomial in latitude and then produced the two-per-second (2-HZ) SSHs in the 2 nd second from the solved polynomial coefficients. Note that the resampled 2-HZ SSHs are located in the 2nd second of the 3-s arc. The use of 3-s data for fitting the 20-HZ SSHs results in significantly smoother 2-HZ SSHs, compared with the use of 1-s data for fitting. Pope's (1976) tau-test procedure was used simultaneously to detect outliers when resampling the 20-HZ data.
In the resampling, each of the 3-s SSHs (60 samples) is expressed as a quadratic function of latitude where l i is the observed SSH, v i is the residual,ϕ i is the latitude of l i , and a, b and c are the coefficients of the quadratic polynomial. Equation (1) can be represented in the following matrix form The least-squares solution of the vector of coefficients where P is the weight matrix, which is an identity matrix in this study. Then the residual vector can be computed as The a posteriori variance of unit weight σ 2 0 is where r n − q is the degree of freedom. For 60 samples and the two-degree polynomial, n 60, q 3. The cofactor matrix of the residuals is An SSH is flagged as an outlier if its residual satisfies the condition (Pope 1976 where q kk is the kth element of the diagonal elements of Q V , 1 − α is the confidence level, τ 1−α; 1, r −1 can be computed by F-distribution F 1−α; 1, r −1 with 1 and r − 1 degrees of freedom. Not all SSHs whose residuals satisfying the condition in Eq. (7) are outliers because large residuals of some SSHs may be caused by a biased solution of X caused by real outliers. Thus, we used an iterative produce to detect the real outliers correctly. In each iteration, we found the maximum of the τ values satisfying the condition in Eq. (7) and removed the corresponding SSH. We repeated this process until we no longer found outliers. After removing all the outliers, the remaining SSHs were approximated by a quadratic polynomial (Eq. 1), and the needed 2-HZ SSHs (at 1.25 and 1.75 s of the 3-s arc) were then computed from the leastsquares-estimated polynomial coefficients in X (Eq. 3) by are the expected latitudes of 2-HZ SSHs, L 2HZ is a vector containing the resampled 2-HZ SSH. According to the error propagation theory, the error variance matrix of the estimated 2-HZ SSHs from the error covariance matrix of the coefficients of the quadratic polynomial X is where X σ 2 0 (A T PA) −1 .

Computing along-track geoid gradients
In this paper, we use the GGs to derive the marine gravity anomaly and VGG by the method of inverse Vening-Meinesz (IVM) and by numerical differentiations, respectively. Because it is not possible to carry out global integrations when converting GGs to gravity anomalies and VGGs, the remove-compute-restore procedure is used (Forsberg 1984). First, the dynamic oceanic topography (DOT) from Levitus et al. (1997) is removed from the altimeter SSHs to generate the geoid heights (GHs). Second, the reference GHs of the Earth Gravitational Field Model 2008 (EGM2008; Pavlis et al. 2012) are also removed to obtain the residual GHs. The standard errors of the EGM2008 coefficients are considered in the global covariance function we used (see Eq. (23) below). Then the residual along-track GGs are approximated by the slope of two successive residual along-track GHs where the ε α, res is the residual GG with azimuth α and has an opposite sign to DOV, N res1 and N res2 are two successive residual GHs, and d is their distance. The full along-track GGs can be obtained when the residual GHs are replaced by the GHs in Eq. (10). The residual GGs are used to compute the north and east components of GGs by LSC described in Sect. 3.2, and the full GGs are used to calibrate the error variances of the GGs from multiple altimeters in Sect. 3.1. The estimated error variance of ε α, res is simply where σ 2 N 1 and σ 2 N 2 are the error variances of N res1 and N res2 , and are the error variances of the corresponding 2-HZ SSHs computed by Eq. (9). The error variance of the full GGs is the same as that of the residual GGs. Because the correlation coefficients of two successive residual GHs (or SSHs) cannot be properly estimated, we simply assume that they are zero to avoid damages introduced by the use of improper correlations, i.e., N res1 and N res2 are uncorrelated. Therefore, in Eq. (11) the covariance of N res1 and N res2 is ignored.

Shipborne gravity anomalies
The shipborne gravity anomalies in the northern part of the SCS are for validating the marine gravity anomalies derived in this paper. From the Ocean Data Bank of Taiwan, we obtained the shipborne gravity anomalies measured by Taiwan's R/V Ocean Researcher I (ORI) over 2009-2016. The standard deviation and the mean value of the crossover differences of the gravity anomalies are 2.2 and 2.3 mgal, respectively (Lee et al. 2016). Figure 1 shows the distribution of the shipborne gravity data. For assessing the gravity anomalies in this paper, we selected high-accuracy shipborne data along 19 straight-trajectory lines (L1-L19 in Fig. 1). The 19 lines are straight so that the influence of the acceleration due to turning can be avoided, and they are located in areas with rough gravity fields, having short-wavelength gravity signals. The fluctuations of the gravity anomalies along the 19 lines are small, suggesting that the gravity data are stable. Lee et al. (2016) used the same shipborne gravity data processing software as the one used in . Thus, the spatial resolution of the along-track shipborne gravity anomalies is about 0.5 km ), corresponding to 0.25', compared with the 1' grid interval for the altimeter-derived gravity anomalies in this paper. The shipborne gravity data is affected by long-wavelength errors such as mechanical drift of the gravimeter and absence of base-station ties (Wessel and Watts 1988). We used the gravity field of EGM2008 (Pavlis et al. 2012) to degree 2160 as a reference field to remove the systematic bias of the shipborne measurements of the 19 lines. For each line, the differences between the shipborne gravity and the reference gravity field were fitted by a quadratic polynomial of time, and the polynomial coefficients were estimated by the least-squares method. Then the systematic biases determined by the estimated coefficients were used to correct the shipborne data (Hwang and Parsons 1995).

Calibrating the error variance components of multi-altimeter geoid gradients by MINQUE
This study uses the method of LSC with noise (Moritz 1980) to estimate the north and east gradient components from the along-track GGs of the three altimeter missions. In this method we take care of two issues: calibrating the error variances of the GGs from multiple altimeter datasets (this Section) and scaling the global covariance functions (see Sect. 3.2 for the detail of scaling). The LSC method for generating the north and east gradient components using along-track GGs has been presented by Hwang and Parsons (1995). This study advances this method by optimally calibrating the data weights and global covariance functions (see Sect. 3.2). A priori (initial) error variances of the GGs can be assigned empirically or propagated from the standard  Taiwan. L1-L19 are 19 straight-trajectory lines for assessing the gravity anomalies errors of the SSHs that are used to generate the corresponding GGs. Such a priori variances may not be correct weights when data from multiple satellites are used. Therefore, we should calibrate the initial variance to obtain appropriate relative weights of the data from these three altimeter missions. For the variance calibration, this study uses an approach based on MINQUE (Rao 1971;Koch 2013), which determines the variance components of GGs from Cryosat-2, Jason-1, and Jason-2. An along-track GG within a data selection window (Fig. 2) forms an observation equation as where ε i is the observed full along-track GG, v i is the residual, α i is the azimuth of ε i , ξ and η are the north and east gradient components at the grid point, respectively, (ϕ i , λ i ) and (ϕ 0 , λ 0 ) are the latitudes and longitudes of ε i and the grid point, respectively, ϕ i ϕ i − ϕ 0 , and λ i λ i − λ 0 . The term ξ cosα i + ηsinα i is the sum of the projections of the north and east gradient components in the along-track direction (Heiskanen and Moritz 1985).
f is used to absorb the spatial variation of the GGs within the data selection window. The window should be big enough to have sufficiently large degrees of freedom for least-squares estimating the parameters in Eq. (12), and it should also be small enough to well satisfy the point relationship between the GGs and their component, i.e., ε i ξ cosα i + ηsinα i .
Both Jason-1 and Jason-2 are the follow-on missions of the TOPEX/POSEIDON mission and the two missions result in SSHs with similar standard errors (Sandwell et al. 2019). In the low-resolution model (LRM), Cryosat-2's altimeter functions like a conventional radar altimeter such as the Jason-1 and Jason-2 altimeters but produces more accurate SSHs than the latter two. Therefore, we divide the altimeter data into two categories, i.e., those from Cryosat-2 and from Jason-1/2. The observation equations in Eq. (12) for the GGs from these two categories (Cryosat-2 and Jason-1/2) in a data selection window like that in Fig. 2 are represented by the following matrix forms Let V V 1 V 2 be the joint vector of residuals, A be the joint design matrix, i.e., A A 1 A 2 , and L L 1 L 2 be the joint vector of observations with L 1 and L 2 containing GGs from Cryosat-2 and Jason-1/2, respectively. Equation (13) can be written as If a data selection window (Fig. 2) The elements in the vector of unknows X are X ξ η a b c d e f T . We assume that these two categories of GGs are independent of each other so that the covariance matrices between vectors L 1 and L 2 are zero matrices. Hence, the initial covariance matrix for all the Fig. 2 Diagram of determining the north and east gradient components (ξ, η) at a grid point (ϕ 0 , λ 0 ) from along-track GGs (colored dotted lines representing different altimeter datasets). The green box shows the data selection window GGs is where C 1 and C 2 are the cofactor matrices of GGs from Cryosat-2 and Jason-1/2, respectively, σ 2 1 and σ 2 2 are the variances of unit weight, which are unknown variance components and will be estimated, and T i is constructed by C i and is an augmented matrix with the same dimension as 0 .
The principle of MINQUE is summarized below. In our case, an arbitrary linear function of the two variance components is where Y contains the two variance components. can be estimated by a quadratic form of the observation vector L, i.e., where M is an unknown symmetrical matrix. In the MINQUE theory, must satisfy the following three properties: (1) Invariance. The quadratic form L T ML is independent of the vector of unknows X, which is equivalent to (2) Unbiasedness. The quadratic form L T ML should be the unbiased estimate of . That is, M should satisfy the following condition.
(3) Minimum norm. The Euclidean norm of the difference between the actual and the theoretical estimated values of Eq. (17) should be minimum, that is, The essence of the MINQUE theory is to estimate the matrix M that minimizes tr(M 0 M 0 ), subject to the constraints in Eqs. (18) and (19). The detailed derivations are given by Rao (1971) and Koch (2013). Then the resulting matrix M is used to estimate the variance components.
By estimating the variance components, we refined the relative weights of the along-track GGs from Cryosat-2 and Jason-1/2. The initial covariance matrix 0 in Eq. (15) can be obtained by setting σ 2 1 σ 2 2 1 with C 1 and C 2 being cofactor matrices obtained from the initial error variances given in Eq. (11). In practice, at each grid point, we set a 8 × 8 data selection window (Fig. 2). Then the variances of the GGs for Cryosat-2 and Jason-1/2 propagated from the error variances of the 2-HZ SSHs form the diagonal elements of the initial cofactor matrices C 1 and C 2 , respectively. A summary of the formulae used in the MINQUE estimation of variance components is as follows (in the 8 × 8 data window around a grid point): The calibrated variance matrix for GGs: 0 σ 2 1 C 1 0 0 σ 2 2 C 2 (6) Replace 0 by 0 and repeat steps (1)-(5) until σ 2 1 σ 2 2 In this study, because we assume that the a priori covariance matrix 0 and the a posterior covariance matrix 0 are both diagonal matrices, they are actually variance matrices. The W matrix in the above procedure is originated from the least-squares estimation of X as Furthermore, the 2 × 2 covariance matrix of the estimated variance components ( σ 2 1 and σ 2 2 ) is where σ 0 is the a priori standard error of unit weight. The diagonal elements of Y are the error variances of σ 2 1 and σ 2 2 . As noted by Caspary (1988, p. 100), a successful estimation of variance components depends on the following two conditions: (1) the data are not colinear and the error effects are separable, (2) the function model (Eq. 12) is suitable so that the groups of data are correlated and there is a sufficiently large data redundancy. Because the inclination angles of Cryosat-2 (92°) and Jason-1/2 (66.04°) are significantly different, Condition 1 is met. Condition 2 is likely to be met because the data densities from these missions are sufficiently large (see Sect. 5.2). The refined error variance matrix 0 is the noise variance matrix in the method of LSC (see Sect. 3.2) for estimating the north and east gradient components.

Scaling the global covariance function
The second issue, entangled with the error variance calibration, is scaling the global covariance functions. First, the global covariance function of the earth's anomalous potential relative to a gravitational model such as EGM2008 can be expressed as (Tscherning and Rapp 1974;Hwang and Parsons 1995) (23) where ψ is the spherical distance, E 2 n is the anomaly error degree variance of degree n of the reference gravity field (from the error variances of the coefficients of EGM2008 gravitational model, Pavlis et al. 2012), R is earth's mean radius (about 6371 km), P n (cosψ) is the Legendre polynomial, and c n A(n−1) (n−2)(n+B) s n+1 is the Model 4 "signal" variance of global gravity anomaly field, with A 425.28 mgal 2 , B 24, and s 0.999617 (Tscherning and Rapp 1974). The first term in Eq. (23) is the error part and the second term is the signal part of K (ψ).
By covariance propagation (Moritz 1980), the global covariance function K (ψ) can be used to construct the covariance functions for the longitudinal and transversal GGs (Tscherning 1974) and then for along-track GGs (Hwang and Parsons 1995, Appendix A). These global covariance functions can be used in the LSC formula to compute the residual north and east gradient components (relative to EGM2008) (Moritz 1980, p. 102): where vector L contains the residual along-track GGs calculated by Eq. (10) within the data selection window in Fig. 2, s is a vector containing the residual north (ξ ) and east (η) gradient components, which are simply the gradients with azimuths of 0°and 90°and are two special cases of along-track GGs, 0 is the calibrated error variance matrix in Sect. 3.1 and works as a filter, C sL is the covariance matrix for signals (in vector s) and observations (in vector L), and C LL is the covariance matrix for observations. The elements in C sL and C LL are constructed from the modified global covariance functions by scaling factors for the residual alongtrack GGs described below. On each grid point, the residual north and east gradient components are computed by LSC using all the residual along-track GGs in the data selection window (Fig. 2). Because we use a small data selection window (8 × 8 ) in this study, we only calibrate the variances of the global covariance functions without calibrating parameters such as correlation length and curvature parameter that characterize local covariance functions. Note that there are several methods for modifying global covariance functions into local covariance functions, e.g., Knudsen (1987), or for directly constructing local covariance functions for a given study area, e.g., Forsberg (1987). As shown in Eq. (23), the global covariance function K (ψ) contains the errors in the coefficients of EGM2008 (error) and the Model 4 degree variances (signal) and is independent of the actual roughness of the gravity field in a data selection window around the grid point (Fig. 2) where the north and east gradient components are estimated. In LSC, the error variances work as a filter, but this filter can function properly only when there is a best balance between signals (for constructing the elements in C sL and C LL ) and data errors (for constructing the elements in 0 ). Otherwise, LSC may result in poor qualities in the estimated north and east gradient components.
In this paper, we use the following approach to modify the global covariance functions into local covariance functions. Specifically, in a given data selection window (Fig. 2), we determine a scaling factor for the global covariance function using: where V 2 ξ and V 2 η are the variances of the residual north and east gradient components, respectively, computed by Eq. (24) in the data selection window, and C ll (0) and C mm (0) are variances of longitudinal and transversal gradients derived from the global covariance function (the spherical distance is zero) in Eq. (23). In this study, the global covariance functions for the longitudinal and transversal gradients are computed by a FORTRAN program coded by Tscherning (1974). The idea of scaling global covariance function has been proposed in previous studies, which were mostly based on sparse altimeter data, e.g., Hwang et al. (1998). The present altimeter data density is vastly increased over that in the 1990s, thus we expect a different result from the use of scaled global covariance functions for optimized marine gravity determinations from altimetry.
To grid the scaling factors (β), we also use a moving window to compute the β values on the central grid. We experimented with three window sizes, i.e., 10 × 10 , 20 × 20 , and 30 × 30 . Using the results of shipborne gravity assessments, we found that the resulting gravity anomalies differ by only 0.001 mgal. One criterion for selecting a proper window is that the computed scaling factors should be sufficiently smooth from one grid point to the next grid points (in the four quadrants). For the best smoothness and smallest computational load, we decided that 20 × 20 is the best window size. We then computed scaling factors on a 1 × 1 grid (the same grid for gravity anomaly and VGG). From visual inspections, the scaling factors obtained in this way can avoid oscillating north and east gradient components in the subsequent LSC computations (see Eq. 26 below).
Within a small spherical distance, the global covariance function K (ψ) and the global covariance function of alongtrack GGs roughly are nearly Gaussian (Hwang and Parsons 1995). Thus, the needed local covariance functions for the matrices in Eq. (24) can be obtained by multiplying the global covariance functions by β. As such, we can re-write Eq. (24) as where C sL and C LL are scaled covariance matrices, which are constructed from the scaled global covariance functions, i.e., β K (ψ), and the covariance functions propagated from it. The simultaneous use of calibrated error variances of GGs and scaled global covariance functions as in Eq. (26) can balance the contributions from signals and noises to improve the estimated north and east gradient components.

Gravity anomaly
The gridded north and east gradient components are used to derive gravity anomalies and VGGs. The methods are summarized in the two subsections here. The regular grids of residual north and east gradient components generated in Sect. 3.2 can be converted to residual gravity anomalies by the inverse Vening-Meinesz (IVM) formula based on 1D FFT method (Hwang 1998): where P is the computational point, Q is a contributing point on the unit sphere, g P is the residual gravity anomaly at P, γ 0 is the normal gravity, ξ cos ξ cosϕ, η cos ηcosϕ, λ Q P λ Q − λ P , ϕ and λ are grid intervals in the directions of latitude and longitude, respectively, α Q P is the is a kernel function depending on spherical distance (ψ Q P ), and F 1 is the 1D FFT operator. By adding the reference gravity values of EGM2008 to the residual gravity anomalies computed by Eq. (27), the full marine gravity anomalies can be obtained.
We also add the innermost-zone effect around the neighborhood of the computational point P, where the kernel function H is nearly or completely singular. The innermostzone effect on the gravity anomaly is where ξ y ∂ξ /∂y and η x ∂η/∂x (x and y are the rectangular coordinates pointing east and north, respectively) are the gradients of the GGs, and x and y are the grid intervals. ξ y and η x are determined by numerical differentiations of the gridded ξ and η values and can be used to determine VGGs (see Sect. 4.2). Figure 3 shows the flowchart of deriving marine gravity anomalies from multi-altimeter data in this paper. Fig. 3 The flowchart of deriving marine gravity anomalies from multi-altimeter data

Vertical gravity gradient
The VGG in this study is the vertical gradient of the gravity anomaly and has the same definition as that given by Heiskanen and Moritz (1985) and Sandwell and Smith (1997, except the sign). In fact, a gravity gradient is the sum of the gradient of the normal gravity defined by an ellipsoid and the gradient of the gravity anomaly. The VGG in this study and in Sandwell and Smith (1997) is only the vertical gradient of the gravity anomaly. We calculate the VGG from the gradient components obtained in Sect. 3.2. The formula for determining VGG has been given by Heiskanen and Moritz (1985) and includes more terms than the formula given by Sandwell and Smith (1997). According to Heiskanen and Moritz (1985, p. 117), the VGG in the spherical coordinate is where g is the gravity anomaly, r is the radial variable, which is positive toward the opposite direction of gravity, G is the mean gravity value over the earth, R 6371 km, is earth's mean radius, N is the geoidal undulation, λ and ϕ are longitude and latitude, respectively, ξ and η are the north and east DOV components and have opposite signs with the north and east gradient components (ξ and η) in Sects. 3.2.
Using the same local rectangular coordinates (x, y) as in Eq. (28) and the relationships Rdϕ dy, Rcosϕdλ dx, we have ∂ g ∂r According to Heiskanen and Moritz (1985), the first two terms on the right-hand side are small compared to the third term; this issue will be discussed in Sect. 5.5. The third term in Eq. (30) can be readily determined from the numerical differentiations of the north and east gradient components for the innermost-zone effect (Eq. 28). The residual north and east gradient components (ξ and η) determined by Eq. (26) are used for the numerical differentiations to generate the term for the innermost-zone effect (Eq. 28) and the third term in Eq. (30). The reference VGG is obtained from the reference field of the DOV components (ξ and η ) based on the EGM2008 model to degree 2160. By restoring the reference VGG to the residual VGG, the full VGG can be obtained. Note that the sign of DOV components of EGM2008 is opposite to that of gradient components defined in Sect. 3.2.

Resampled 2-HZ data
Using the procedure described in Sect. 2.1.2, we resampled the 20-HZ SSHs into 2-HZ SSHs. As an example, Fig. 4 shows the result of resampling 3-s SSHs along a section of pass 127, cycle 501 of Jason-2. In this example, only one outlier was detected and removed, and the remaining data were fitted by the quadratic polynomial. Then the 2-HZ SSHs were obtained at the 2nd second along the fitted polynomial (the green line, circles are the 2-HZ SSHs at 1.25th and 1.75th second). The fitted quadratic polynomial from the 3-s data with/without the outliers is shown by the cyan/black line. The outlier has affected the resampled 2-HZ SSHs but not to a serious extent. In the study area shown in Fig. 1, there are 6,935,601 20-HZ Cryosat-2 SSHs, and a total of 258,837 outliers (3.73%) were detected. For Jason-1/2, the numbers of outliers are 31,798 (5.73%) and 54,779 (5.84%), respectively. The outliers were not used in the determination of geoid gradients. The shortest wavelengths corresponding to the resampled 2-HZ SSHs of Jason-1/2 and Cryosat-2 are 5.8 km and 6.4 km, respectively. The mean cross-track distances of Jason-1 GM, Jason-2 GM and Cryosat-2 at the equator are 8 km, 4 km and 7.5 km, respectively.
It is possible that the 2-HZ SSHs resampled from 3-s data may be over smoothed, resulting in the loss of highfrequency seafloor characteristics. To examine this issue, we also generated 2-HZ SSHs by fitting the 1-s (20 samples) data with a quadratic polynomial. Figure 5a shows the 2-HZ SSHs generated by fitting 3-s (60 samples) and 1-s (20 samples) data of Jason-2. The 2-HZ SSHs generated from the 3-s data are smoother than those from the 1-s data. The along-track gradients of the 2-HZ SSHs from the 3-s data are smoother than those from the 1-s data, as shown in Fig. 5b. Note that smooth SSHs and gradients do not necessarily lead to more accurate gravity anomalies and VGGs. To show the effectiveness of resampled 2-HZ SSHs from the 3-s fitting, we determined the root-mean-squared (RMS) differences between the shipborne gravity anomalies and the Jason-2-derived gravity anomalies by using the 2-HZ SSHs from the 1-s and 3-s fittings. The differences are shown in Table 2. An improvement rate in Table 2 is calculated by rate × 100%, and a positive rate value suggests improvement. Except line 5, the accuracies along all lines are improved by using the 3-s fitting. The deterioration in the accuracy of the 3-s fitting along line 5 may be caused by systematic biases in the 3-s case in some sections of this line. The average RMS differences in the 1-s and 3-s fittings are 4.05 mgal and 3.39 mgal, respectively. Thus, the use of the 3-s fitting improves the gravity accuracy by 15.3% over the use of the 1-s fitting. These results indicate that the resampled 2-HZ SSHs from the 3-s data outperform the case of 1-s data, resulting in marine gravity anomalies with higher accuracies. Because the 2-HZ SSHs from the 3-s data are already smooth, no filtering is applied to such 2-HZ SSHs in order to retain the high-wavenumber seafloor signatures in the altimeter-derived gravity and VGG fields.

Calibrated error variances by MINQUE
The initial error variances of residual along-track GGs for Cryosat-2 and Jason-1/2 are shown in Fig. 6a, c, respectively. They are computed by Eq. (11) using the error variances of the resampled 2-HZ SSHs by Eq. (9). The initial gradient variances of Jason-1/2 are larger than those of Cryosat-2 because Cryosat-2 has an overall higher observation accuracy than Jason-1/2 (Sandwell et al. 2019). In theory, the error variances of GGs should be larger in the coastal regions, where radar waveforms can be contaminated (Deng et al. 2002) and geophysical correction models for SSH measurements can  (RMSD, in mgal) between the shipborne and the Jason-2-derived gravity anomalies (1-s and 3-s fittings) and improvement rate from 1-ec to 3-s fitting be poor. However, Fig. 6a, c show that this is not the case for both categories of GGs. This may be because the original 20-HZ SSHs have been free from waveform contaminations and poor geophysical corrections, the variances of the resampled 2-HZ SSHs are already uniform over coastal regions and open oceans. However, the 20-HZ SSHs in the coastal regions may contain long-wavelength systematic errors, but they do not propagate into the errors of GGs because of the use of along-track differentiations when deriving GGs. As such, the variances of GGs in Fig. 6a, c are mostly caused by the random errors of the 20-HZ SSHs not by the potential systematic errors of the SSHs, especially in the coastal regions. Using the MINQUE theory described in Sect. 3.1, we estimated new error variances of along-track GGs from Cryosat-2 and Jason-1/2 for each grid point in a 8 × 8 data selection window. The calibrated (new) error variances are shown in Fig. 6b, d. The patterns of the calibrated variances (Fig. 6b, d) are different from those of the initial variances (Fig. 6a, c). The calibrated variances tend to increase with the roughness of the marine gravity field. For example, the calibrated variances are relatively large around the Dongsha Island and the area to its east, the Batan Islands, and the Babuyan Islands and their coastal areas. One reason for this relationship is that the calibration by MINQUE will be Fig. 6 Left column: the initial error variances of residual GGs for a Cryosat-2 and c Jason-1/2, and e their ratios (variances of Jason-1/2 divided by variances of Cryosat-2). Right column: the calibrated error variances of residual GGs for b Cryosat-2 and d Jason-1/2 by MINQUE, and f their ratios effective only when the along-track GGs are with sufficiently large magnitudes; this occurs only over an area with a rough gravity field where the GGs are relatively large for effectively determining the variance components. The results in Fig. 6b, d suggest that the calibrated variances of Jason-1/2 are larger than those of Cryosat-2. The calibrated error variances are used for an optimal LSC determination of the north and east gradient components by Eq. (26). In addition, the along-track patterns in Fig. 6a, c, e are caused by large error variances of data along several tracks. The MINQUE calibration uses data in a 8 × 8 window, thus reducing such track patterns (Fig. 6b, d, f). Figure 6e shows the ratios between the initial variances of Cryosat-2 and Jason-1/2 (variances of Jason-1/2 divided by variances of Cryosat-2) within the 8 × 8 data selection windows. The ratios are proportional to the numbers of GG observations from Jason-1/2 (Fig. 7d). Figure 6f shows the ratios of the calibrated variances of the two categories of altimeter data, which are smaller than those in Fig. 6e. This is because the calibrated variances of Jason-1/2 Fig. 7 The accuracies (variances, in sec 4 ) of estimated variances components for a Cryosat-2 and b Jason-1/2; the numbers of gradient observations in a 8 ×8 data selection window for c Cryosat-2, d Jason-1/2, and e combined Cryosat-2 and Jason-1/2; and f the distribution of the numbers of iterations for estimating variance components using MINQUE (Fig. 6d) in some regions are smaller than the initial variances (Fig. 6c). This result indicates that the relative weights of the 2-HZ SSHs of Jason-1/2 in some regions are improved after calibrating by MINQUE. The comparison between the sub-figures in the left column of Fig. 6 and those in the right column suggests that using the initial variances of the two categories of altimeter data is inappropriate, and the initial variances of Jason-1/2 are artificially too big in comparison to those of Cryosat-2. Figure 7a, b shows the accuracies (variances) of estimated variance components (Eq. 22) for Cryosat-2 and Jason-1/2, respectively, which are proportional to the numbers of points in the data selection windows (Fig. 7c, d). As such, the accuracies are low in the areas with sparse altimeter data points. In most windows, the numbers of Cryosat-2 are larger Fig. 8 Distribution of the scaling factors of the global covariance functions. Relatively large scaling factors occur in areas of rough gravity fields and shallow waters than 120 (Fig. 7c) and those of Jason-1/2 are larger than 40 (Fig. 7d). Thus, the total number of Cryosat-2 and Jason-1/2 can be larger than 150 (Fig. 7e). The degrees of freedom are sufficiently large to estimate reliable variance components (Fig. 6b, d). Moreover, 2-6 iterations are needed to reach convergences in the variance component estimations, as shown in Fig. 7f.

Scaling factors for the global covariance functions
Using Eq. (25), we computed the scaling factor (β) of global covariance functions from the north and east gradient components estimated by Eq. (24) using the calibrated error variances of GG. Figure 8 shows the scaling factors, which are for modifying the global covariance functions into local covariance functions to improve the estimated north and east gradient components by LSC. The distribution pattern of the scaling factors (Fig. 8) is correlated with the marine gravity pattern (see Sect. 5.4). In general, scaling factors are large in areas where the gravity anomaly variances are large. This is because a rough gravity field corresponds to a large variance of GGs, leading to a large scaling factor that in turn increases the global covariance functions. The global covariance functions are based on the standard errors of the EGM2008 geopotential coefficients and the model gravity signal components beyond degree 2160 (Eq. 23), thus such functions cannot really reflect the local gravity field variations, especially in areas with seamounts, trenches, and other high-wavenumber seafloor signatures. On the other hand, a scaling factor less than one leads to smaller scaled global covariance functions, resulting in an artificial increase in the error variances of GGs in a data selection window and leading to relatively smooth north and east gradient components from LSC. The use of the scaling factors has improved the gravity accuracy (see Sect. 5.4).

Improved marine gravity accuracy due to calibrated error variances and scaled global covariance functions
Using both the calibrated error variances (Fig. 6c, d) and the scaling factors (Fig. 8), we derived marine gravity anomalies on a 1 ×1 grid in the SCS from Cryosat-2 and Jason-1/2 data by IVM, as shown in Fig. 9a. We assessed the accuracy of the altimeter-derived gravity anomalies using the shipborne gravity anomalies along the 19 lines in Fig. 1. The RMS differences are shown in Table 3. The rate in Table 3 is defined as rate RMSD initial −RMSD calibratedandscaled RMSD initial × 100%, which indicates the improvement rate due to calibration and scaling. On average, the use of the calibrated error variances and scaled covariances improves the gravity accuracy by 9.5%, compared with the case using the initial error variances. For line 10, the negative rate may be due to systematic biases in some sections of this line between the gravity anomalies from the altimeter case using the calibrated and scaled covariance and the shipborne measurements. Table 3 also shows that the gravity accuracy from the error-calibrated Cryosat-2 and Jason-1/2 altimeter data and the scaled covariances is better than the gravity accuracy from only un-calibrated Cryosat-2 or Jason-1/2 altimeter data. Thus, the gravity accuracy benefits from the use of multiple altimeters and from error calibration and covariance scaling. The differences between the altimeter-derived gravity anomalies (Fig. 9a) and the shipborne gravity anomalies along the 19 lines in Fig. 1 are shown in Fig. 9b. The differences (Fig. 9b) appear to have no specific patterns. Fig. 9 a Gravity anomalies derived from Cryosat-2 and Jason-1/2 altimeters using the calibrated error variances and scaled global covariance functions, and b the differences between the gravity anomalies in (a) and the shipborne gravity anomalies along the 19 lines in Fig. 1

Vertical gravity gradient determination
VGGs on a 1 × 1 grid were determined by Eq. (30) using the north and east gradient components estimated from the calibrated variances and scaled global covariance functions. Heiskanen and Moritz (1985) recommended that the first two terms in Eq. (30) can be neglected because they are relatively small compared with the third term. However, in this paper, we also computed the first two terms in an attempt to determine optimal VGGs. First, we validate our computer program for VGG computation by comparing the VGGs from this paper and from Sandwell et al. (2014) in the Gulf of Mexico, where extinct ridges have been discovered in the VGGs from Cryosat-2 and Jason-1 altimeter data. Note that, to our knowledge, there are no publicly available observed VGGs in the world oceans for assessing the accuracy of altimeter-derived VGGs. The altimeter-derived VGGs may be filtered and compared with the vertical gradients from the satellite mission GOCE such as those presented by Ebbing et al. (2018). However, this comparison is not carried out in this study and can be a subject of future studies. Figure 10a shows the VGGs in the Gulf of Mexico from the first term in Eq. (30), i.e., 2G R 2 N , where N is the geoidal height of EGM2008 to degree 2160. Figure 10b, c show the reference and residual VGGs from the second term in Eq. (30), i.e., G R ξ tanϕ or − G R ξ tanϕ, respectively. The reference VGGs are from the north DOV components (ξ ) of EGM2008 to degree 2160, and the residual VGGs are from the north gradient components (ξ ) computed in this study. The reference and residual VGGs from the third term in Table 3 RMS differences (in mgal) between the shipborne gravity anomalies and the gravity anomalies derived from Cryosat-2 and Jason-1/2 using the initial variances and the calibrated and scaled covariances

Line
The initial variances The calibrated and scaled covariances Eq. (30), i.e., −G ∂ξ ∂ x + ∂η ∂ y or G ∂ξ ∂ x + ∂η ∂ y , are shown in Fig. 10d, e, respectively. The VGGs in Fig. 10d are computed from the DOV components of EGM2008, and those in Fig. 10e are computed from the residual GGs in this paper. Although the VGGs from the first two terms (Fig. 10a-c) are small in comparison to those from the third term (Fig. 10d, e), they are correlated with seafloor signatures in the Gulf of Mexico. Figure 10f shows the full VGGs generated by combining the quantities in Fig. 10a-e. The tectonic signatures that have been interpreted by Sandwell et al. (2014), such as the extinct spreading ridges, fracture zones, and continent-ocean boundaries, are evident in Fig. 10f and are similar to those indicated in the VGG field of Sandwell V31.1 (Fig. 10g) (the method is documented by Sandwell and Smith (1997); see also https:// topex.ucsd.edu/pub/global_grav_1min/). This shows that our computer program correctly determines the VGGs. Note that the tectonic signatures are mainly contributed by the reference VGGs from the third term (Fig. 10d), with Fig. 10b also showing some tectonic signatures. The residual VGGs from the second and third terms (Fig. 10c, e) show signatures that are due to altimeter data noises and high-wavenumber seafloor variations. Note that the sign of the VGG derived by our program is opposite to that of Sandwell V31.1 VGG.
For a convenient comparison, we reverse the color bars in the VGGs in Fig. 10a-f. In this paper, the sign convention of VGG follows that used in physical geodesy (Heiskanen and Moritz 1985), thus GG has the opposite sign with DOV. Sandwell and Smith (1997) may have treated along-track GGs as along-track DOV. As such, our VGGs have an opposite sign with the VGGs defined in Sandwell V31.1. The VGG from this study (Fig. 10f) and from Sandwell V31.1 (Fig. 10g) have almost identical signatures.
After the validation in the Gulf of Mexico, then we computed VGGs in the northern part of the SCS (same coverage as in Fig. 9a). Figure 11a shows the full VGGs in the northern part of the SCS from the calibrated Cryosat-2 and Jason-1/2 altimeter data using our computer program. The tectonic signatures in Fig. 11a are similar to those found in the VGG field of Sandwell V31.1 (Fig. 11b). Note that some stripes are present in Fig. 11b (Sandwell V31.1), but not in Fig. 11a (this paper). Both the gravity field (Fig. 9a) and the VGG field (Fig. 11a) are correlated with seafloor signatures in the same way, but the VGG field shows sharper boundaries of features than the gravity field. The differences between the VGGs from this paper and from Sandwell V31.1 are shown in Fig. 11c, which appear to be random. The mean and RMS Fig. 10 VGGs in the Gulf of Mexico according to Eq. (30) from a the first term, b the reference and c residual parts of the second term, d the reference and e residual parts of the third term; f the full VGGs, and g Sandwell V31.1 VGGs. Note that the color bar in (g) is opposite to that in (a)-(f) Fig. 11 The full VGGs in the northern part of the SCS from a this paper, b Sandwell V31.1, and c the differences between (a) and (b). Note that the two color bars have opposite color scales because of the opposite signs of VGGs in (a) and (b). The red line at 19.83°N shows a profile along which gravity anomalies, VGGs, and depths are compared in Fig. 12  Fig. 12 a The gravity anomalies (Fig. 9a) and VGGs from this paper (red line, Fig. 11a) and from Sandwell V31.1 (cyan line, Fig. 11b), and b the multi-beam depths along the profile in Fig. 11a values of the differences in Fig. 11c are 0.044 Eötvös and 10.160 Eötvös, respectively. Figure 12a compares the altimeter-derived gravity anomalies (Fig. 9a) and VGGs from this study (Fig. 11a) and from Sandwell V31.1 (Fig. 11b) along the profile indicated in Fig. 12a. Figure 12b shows the oceanic depths from multibeam measurements (Yu et al. 2021). There are notable correlations between the gravity anomalies, VGGs, and the depths near the areas marked by the vertical green dashed lines. The four magenta lines mark the locations where the VGGs appear to be more correlated with small seafloor features than gravity anomalies. This higher sensitivity of VGG shows that VGGs have the potential to recover bathymetry at short spatial wavelengths. A geophysical interpretation of the altimeter-derived VGGs in the northern part of the SCS is not pursued here but can be a subject of future studies.

Conclusions
In this study, we obtained 2-HZ SSHs from 20-HZ SSHs of Cryosat-2 and Jason-1/2 and determined along-track GGs and their initial error variances. The initial variances were then calibrated by the method of MINQUE. For each grid point, a scaling factor was determined to scale the global covariance functions. The calibrated and scaled covariances were used to estimate the north and east gradient components, which are then used to derive gravity anomalies and VGGs. In the northern part of the SCS, the accuracy of the gravity anomalies from the calibrated and scaled covariances are improved compared with the case using the initial error variances of GGs. Our method for computing VGGs considers more terms compared to the methods used in previous studies. We validated the computer program for VGG computation in the Gulf of Mexico and then derived VGGs in the northern part of the SCS. Compared to gravity anomalies, VGGs are more sensitive to seafloor signatures at shorter spatial wavelengths and might be more powerful in detecting structural boundaries.
In this study, we only used Cryosat-2 and Jason-1/2 data to calibrate error variances and scale global covariance functions by the methods given in Sect. 3. We expect that the methods are applied to more altimeter datasets to refine the marine gravity anomalies and VGGs presented in Sect. 5. For example, the range measurements by the SARAL/AltiKa altimeter are the most accurate among all ongoing and past altimeters, its SSH dataset can be useful for error variance calibration and gravity accuracy improvement. The future SWOT altimeter mission will provide wide-swath SSH observations with unprecedented spatial resolutions and range accuracies. Calibrating the data errors of the existing missions with respect to the data errors of the SWOT mission is highly important when merging all these altimeter data.