Comparison of Mean Dynamic Topography Modeling from Multivariate Objective Analysis and Rigorous Least Squares Method

Abstract

Filtering methods are usually used to combine the mean sea surface (MSS) and geoid (computable by global geopotential model (GGM)) into a common subspace, to model mean dynamic topography (MDT), which may lead to signal leakage and distortion problems. The use of the rigorous least squares (LS) method and multivariate objective analysis (MOA) alleviates these problems, and the derived MDTs from these two methods show better performance than MDTs derived from filtering methods. However, the advantages and disadvantages of these two methods have not been evaluated, and no direct comparison has yet been conducted between these two approaches regarding the performances in MDT recovery. In this study, we compare the performances of the MOA method with the LS method, providing information with respect to the usability of different methods in MDT modeling over regions with heterogeneous ocean states and hydrological conditions. We combined a recently published mean sea surface called DTU21MSS, and a satellite-only GGM named GO_CONS_GCF_2_DIR_R6, for MDT computation over four typical study areas. The results showed that the MDTs derived from the LS method outperformed the MOA method, especially over coastal regions and ocean current areas. The root mean square (RMS) of the discrepancies between the LS-derived MDT and the ocean reanalysis data was lower than the RMS of the discrepancies computed from the MOA method, by a magnitude of 1–2 cm. The formal error of the MDT estimated by the LS method was more reasonable than that derived from the MOA method. Moreover, the geostrophic velocities calculated by the LS-derived MDT were more consistent with buoy data than those calculated by the MOA-derived solution, by a magnitude of approximately 1 cm/s. The reason can be attributed to the fact that the LS method forms the design matrix segmentally, based on the error characteristics of the GGM, and suppresses high-frequency noise by applying constraints in different frequency bands, which improves the quality of the computed MDT. Our studies highlight the superiority of the LS-derived method versus the MOA method in MDT modeling.

1. Introduction

Mean dynamic topography (MDT) is an important data source for studying land/sea datum unification, sea level change and climate change [1,2,3]. The geostrophic current velocity can be easily calculated by using the derivative of the MDT. The geostrophic current velocity has an important influence on ocean dynamic processes such as material migration and heat transfer and exchange, as well as on human activities [4]. Accurate modelling of ocean currents has important implications for meteorology, oceanography, and geophysics [5,6].

The MDT can be determined by combining a precise geoid (which can be calculated by a global geopotential model (GGM)) and mean sea surface (MSS). At present, the accuracy of the MDT, estimated by combining an MSS and a satellite-only GGM, reaches decimeter level in coastal areas and centimeter level in open sea areas [7]. Due to the different spatial resolution between the geoid (10–100 km) and MSS (a few kilometers), the direct combination of these two data sets results in spectral aliasing and spectral leakage [8]. Thus, the spectrum between the MSS and the geoid needs to be homogenized with a filter. Several filtering methods, such as Gaussian filtering, wavelet filtering, and adaptive filtering, can be used to combine the MSS and the geoid to compute the MDT [9,10,11]. However, these filtering methods still suffer from signal leakage and distortion problems, and the formal error of the associated MDTs cannot be estimated from these methods [12].

To mitigate these problems, other methods have been proposed, such as the rigorous least squares (LS) method [13,14,15] and multivariate objective analysis (MOA) [7,16,17,18]. Becker et al. [13] proposed the LS method to estimate MDT, in which the Lagrange basis function (LBF) was used to parameterize MDT. In the LS method, the MDT can be derived by combining the MSS and geoid in a spectrally consistent way based on the LS system, and the design matrix is constructed piecewise, which introduces the error information of the GGM in different bands [19,20]. Previous studies have shown that the MDT modeled from the LS method was of better quality than the MDT modeled from the filtering method (e.g., the Gaussian filtering), by a magnitude of several centimeters [15]. Moreover, Rio et al. [16,17,18] introduced a multivariate objective analysis (MOA) method to estimate MDT, which considered the covariance between observations and the error of observations. The MDT of a grid point is calculated based on the weighted average of the surrounding observations; the weight is related to the variance and covariance of the observations. The improvement of MDT estimated by MOA method is mainly in short-scale signals, and previous studies have shown that the MOA-derived MDT outperformed the traditional Gaussian-filtering-derived MDT, especially in coastal areas [21]. In addition, the MOA method can combine the raw MDT (MSS minus geoid directly) with external data related to MDT, such as buoy data or ocean model data, to model more detailed signals of MDT.

As mentioned above, both LS and MOA methods can be used to obtain MDT, which outperform the MDT estimated by Gaussian filtering. The LS method and MOA method each have their own advantages in terms of modeling MDT. However, the existing research has not compared the advantages and disadvantages of these two methods, and lacks a direct comparative analysis of these two methods. Moreover, the accuracy of the two methods has not been verified under different marine hydrological conditions. This study focuses on comparing MDT modeling based on the MOA method and the LS method. In particular, we evaluate the performances of these two methods over different oceanic areas with heterogeneous ocean state and hydrological conditions, which can provide the proper choice of modeling approach in computing MDT by merging heterogeneous data sets. The structure of this study is as follows. The principles of the rigorous LS method and MOA method are reviewed in Section 2. In Section 3, four study areas are described, and the datasets for local MDT recovery and validation are also introduced. Then, the numerical experiments are displayed in Section 4. A discussion evaluating the MDT obtained by different methods by exciting MDT and ocean models is shown in Section 5. In Section 6, the conclusions are summarized.

2. Method

2.1. Rigorous Least Squares Method

The rigorous least squares method obtains MDT from MSS and GGM by using an LS-based method. It is crucial to establish a complete observation equation and an associated weight matrix in the LS system, which have an important impact on estimating MDT in the LS method. The sum of the geoid derived from a GGM and MDT is MSS.

$$MSS(\mathsf{\theta },\lambda )=Geoid(\mathsf{\theta },\lambda )+MDT(\mathsf{\theta },\lambda )$$

(1)

where $\mathsf{\theta }$ is the latitude in the spherical coordinate system and $\lambda$ is the longitude. The geoid and MSS models used for MDT modeling should be unified in the same resolution grid. The MSS and the geoid are unified to the GRS80 ellipsoid and tide free system.

The Lagrange basis functions (LBFs) can be applied to parameterize the MDT [20].

$$MDT\left(\overline{\theta },\overline{\lambda }\right)={\displaystyle \sum _{k\in K}{a}_k{b}_k\left(\overline{\theta },\overline{\lambda }\right)}$$

(2)

where $b_k$ represents the basis function; $K$ is the number of basis functions; $a_k$ is the MDT value at $\left(\overline{\theta },\overline{\lambda }\right)$; and $\overline{\theta }$ and $\overline{\lambda }$ represent longitude and latitude at nodes, respectively.

The choice of LBF is important for MDT computation. In this paper, a basis function with 16 parameters (16P) is introduced to parameterize the MDT [16,18]. To reduce the correlation between grids, the grid resolution of MSS and GGM used in this paper is set as 0.5°. In order to calculate an unknown MDT point, LBF interpolation with 16 parameters is performed with 4 surrounding points. Then, all grid points of MDT are parameterized by LBF, and the parameter coefficient matrix (A_mdt) is obtained.

In this study, a satellite-only GGM is used for modeling MDT, its maximum expansion degree and order (d/o) is 300. The GGM expression is divided into three parts for processing, according to the signal to noise ratio (SNR) of the GGM. The SNR of the GGM decreases when the d/o of GGM increases. The first part ($cs1$) represents the GGM signals with high SNR, which has spherical harmonics (SHs) from d/o 2 to a suitable cut-off d/o of GGM (e.g., where SNR > 1). The second part ($cs2$) can be recognized as a buffer between cs1 and cs3. The third part ($cs3$) represents the geoid signals that cannot be obtained in the satellite-only GGM, which has SHs from max d/o of $cs2$ to infinity. More detailed information about the MDT parameterization can be found in [20]. Equation (1) can be expressed as:

$$MSS+v=[J_{cs1} J_{cs2} J_{cs3} J_{mdt}]\left[\begin{array}{l}{X}_{cs1}\\ X_{cs2}\\ X_{cs3}\\ X_{mdt}\end{array}\right]$$

(3)

where $v$ is the residual in the LS system, and $X_{cs}$ and $X_{mdt}$ represent the unknown SH coefficients and the MDT values based on LS theory.

The third part ($cs3$) is usually ignored or set to zero [15], due to the limited spatial resolution of the GGM which lacks the information of $cs3$. The signal of $cs3$ is $S=J_{cs3}\cdot X_{cs3}$. We then set $MSS=MSS-S$, as:

$$MSS+v=[J_{cs1} J_{cs2} J_{mdt}]\left[\begin{array}{l}{X}_{cs1}\\ X_{cs2}\\ X_{mdt}\end{array}\right]$$

(4)

To obtain a slightly smooth MDT, additional smoothing information should be added to the observation equation. The smoothing information that causes the norm of the MDT gradient decreases can be added in the observation equation by:

$$\left[\begin{array}{l}0\\ 0\end{array}\right]+\left[\begin{array}{l}{v}_{md{t}_x}\\ v_{md{t}_y}\end{array}\right]=\left[\begin{array}{l}0 0 \nabla {\mathrm{J}}_x\\ 0 0 \nabla {\mathrm{J}}_y\end{array}\right]\left[\begin{array}{l}{X}_{cs1}\\ X_{cs2}\\ X_{mdt}\end{array}\right]$$

(5)

where $\nabla {\mathrm{J}}_x$ is the derivative of the parameterized MDT in zonal, and $\nabla {\mathrm{A}}_y$ is that in meridian.

The complete observation equation can be expressed as:

$$\left[\begin{array}{l}MSS\\ GG{M}_{cs1}\\ 0\\ 0\\ 0\end{array}\right]=\left[\begin{array}{ccc}{J}_{cs1}& J_{cs2}& J_{mdt}\\ I& 0& 0\\ 0& I& 0\\ 0& 0& \nabla {\mathrm{J}}_x\\ 0& 0& \nabla {\mathrm{J}}_y\end{array}\right]\left[\begin{array}{l}{X}_{cs1}\\ X_{cs2}\\ X_{mdt}\end{array}\right]$$

(6)

$$P=\left[\begin{array}{lllll}{K}_{MSS}^{-1}& & & & \\ & K_{cs1}^{-1}& & & \\ & & K_{cs2}^{-1}& & \\ & & & K_I^{-1}& \\ & & & & K_I^{-1}\end{array}\right]$$

(7)

where $P$ is the weight matrix, and $K$ is the variance information which can be obtained by Kaula’s rule.

Assuming that the error of the observation in the observation equation is Gaussian distribution, the observation equation can be expressed as:

$$L-v=A\cdot x,\hspace{1em}\hspace{1em}E\left\{v\right\}=\mathit{0},\hspace{1em}\hspace{1em}D\left\{v\right\}={\sigma }^{2}Q={\sigma }^2{P}^{-1}$$

(8)

where $L=MDT(\theta ,\lambda )$, $A$ and $x$ represent the coefficient matrix of the basis function and the MDT results, respectively; $E$ and $D$ represents the expectation and variance of the observation equation; $Q$ and $P$ represent cofactor matrix and weight matrix; and ${\sigma }^2$ is variance of unit weight.

2.2. Multivariate Objective Analysis

Bretherton et al. were the first to use the multivariate objective analysis method to recover MDT models [22]. Then Rio et al. developed MOA method [16]. The MOA method can be explained as a weighted average method, whose weight is related to the variance and covariance of observations, which allows the method to retain more detailed signals. The detailed information about the MOA method has been discussed by Rio et al. [16] and Wu et al. [7,21]. Moreover, MDT modeled by MOA method preformed the estimated MDT based on the traditional filter method (e.g., Gaussian filtering) [7].

The MDT estimated using the MOA method is given by:

$$\langle h\rangle (r)={\displaystyle \sum _{i=1}^N{\alpha }_{i}O(r_i),{\alpha }_i={\displaystyle \sum _{j=1}^N{A}_{i,j}^{-1}{C}_{r,j}}}$$

(9)

where $\langle h\rangle$ represents the estimated MDT; $\mathit{r}$ represents the grid point; $O(r_i)$ represents the raw MDT observation that was computed by removing the geoid/quasi-geoid directly from MSS; $A$ represents the covariance matrix of the observations; and $C$ is the covariance vector between the observed and estimated MDT. Under the assumption that MDT is isotropic and homogeneous, the covariance only depends on the distance between the observations and the error of the observations [21].

$$\{\begin{array}{l}A={(\langle {\sigma }^2\rangle C(d_{ij})+\langle {\epsilon }_i{\epsilon }_j\rangle )}_{i,j=1,N}\\ C_r={(\langle {\sigma }^2\rangle C(d_{rj}))}_{j=1,N}\end{array}$$

(10)

where ${\sigma }^2$ is the prior variance of the MDT; ${\epsilon }_i$ represents the MDT prior error at grid point i; and $C(r)$ represents the prior covariance function of MDT. This study adopts the prior covariance function introduced by Arhan and De Verdiére [23].

In order to satisfy the application condition of this method, that the mean of estimated MDT must be zero, the residuals that are obtained by deducting a large-scale MDT from the raw MDT are taken as the observations. The large-scale MDT is calculated by filtering the raw MDT by using a Gaussian filter.

The key parameters of the MOA method are the error of the observations, and the variance and covariance between the observations that are related to the distance of the observations. The brief process of estimating MDT by MOA method is as follows. A large-scale MDT is first calculated by applying a Gaussian filter to raw MDT. Then the residual MDT that is obtained by subtracting the large-scale MDT from the raw MDT is set as the observation of the MOA method. The weight of the MOA method is estimated by the exciting reference MDT (for detailed information, refer to Wu et al. [7]). After the residual MDT is improved by the MOA method, the final reconstructed MDT is obtained by adding the large-scale MDT to the improved residual MDT.

Moreover, the error of the estimated MDT could be obtained by [16]:

$$\epsilon (r)={\sigma }^2-{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{A}_{i,j}^{-1}{C}_{r,i}{C}_{r,j}}}$$

(11)

3. Data and Study Area

Four regions were chosen as study areas: Kuroshio Current (5–45°N, 110–150°E), Mexico Gulf (15–55°N, 40–80°W), Agulhas Current (10–50°S, 10–50°E) and East Greenland Current (30–70°N, 25–65°W), respectively. The information of the study areas is shown in Figure 1 (available at www.shadedrelief.com (accessed on 15 May 2022)). There are ocean currents in these study areas, such as Kuroshio Current, Gulf Current, Agulhas Current and East Greenland Current, which are crucial in global water circulation and climate regulation. For example, the Kuroshio Current carries warm water from the equatorial Pacific to near Japan, which warms the coastal areas of southern and southeastern Japan. [24]. The Gulf Current supports the major fisheries in the United States, Mexico and Cuba. The Gulf Current and the Kuroshio Current have a heavy influence on weather conditions in the northern hemisphere [25]. The Agulhas Current, the largest western boundary current in the world’s oceans, is located in the southwestern Indian Ocean, which is crucial for heat transfer and exchange in the South Atlantic [26]. The Greenland Current area contains several currents, namely, the East Greenland Current, the West Greenland Current, the Labrador Current and the North Atlantic Current. These currents are cold and of low-salinity, and are crucial to the transfer of heat through the Arctic and Atlantic Oceans [27]. To evaluate the performances of the MDTs modeled from different methods, we investigated the characteristics of the MDTs’ profiles (red dashed lines in Figure 1) in three areas with different ocean state and hydrological conditions, i.e., coastal area (profile 1), open sea area (profile 2) and ocean current area (profile 3). Section 3.1, Section 3.2, Section 3.3 and Section 3.4 introduce the data we used.

3.1. Mean Sea Surface Model

The DTU21MSS is the mean sea surface model applied to MDT recovery, which is a newly released model by the Technical University of Denmark (DTU). The spatial resolution of DTU21MSS is 1′ × 1′, and the reference time is from 1993 to 2012. The accuracy of DTU21MSS is about 5 cm in open sea areas. For the derivation of DTU21MSS, multi altimetry satellite data are used, such as T/P, Jason1/2, ERS1/2, Sentinel-3A and Cryosat-2. Compared with the previous generation of MSS, updated altimetry data (such as Sentinel-3A/3B) and a modified waveform retracker has been applied in DTU21MSS [28,29]. These improvements make DTU21MSS more suitable for MDT modeling than previous MSS [7].

3.2. Global Geopotential Model

The geoid we used to model MDT was derived from the satellite-only global geopotential model (GGM), instead of a combined model, because the latter increased the computational workload in using the LS method, which was impossible to compute without a large supercomputer. Therefore, GO_CONS_GCF_2_DIR_R6 (DIRR6) was selected as the GGM in MDT recovery. Compared with the previous generation of GGM, recalibrated GOCE gravity gradients and reprocessed orbits that reduce dynamic orbits are used to calculate the DIRR6, which improves the accuracy of DIRR6. The degree and order (d/o) of DIRR6 can be up to 300, corresponding to a spatial resolution of about 66 km [30]. The error degree variance of DIRR6 in geoid height is below 2 cm within d/o 230. The DIRR6 is more accurate than the previous generation of GGMs.

3.3. Synthetic/Ocean MDT Models

To evaluate the MDT estimated by different methods, we introduced several synthetic/ocean numerical models as reference models, such as Simple Ocean Data Assimilation 3 (SODA3) [31], Ocean Reanalysis System 5 (ORAS5) [32], Copernicus, and CNES-CLS18MDT [33]. SODA3 was established by ocean reanalysis method, which improves model resolution, observations, and forced data. This model collected monthly average ocean data from 1980 to 2017, with a horizontal resolution of 1/4°. ORAS5 is an ocean reanalysis model released by ECMWF, which applies the same a priori ocean models and data assimilation approaches as Ocean ReAnalysis Pilot 5 [32]. The ocean data of ORAS5 is a monthly average spanning from 1979 to 2018 with a horizontal resolution of 0.25°. Copernicus was estimated using the DUACS processing system, providing daily data from 1993 to 2018 with a horizontal resolution of 0.25°. The CNES-CLS18MDT is a new mean dynamic topography model that was released by Centre National d’Etudes Spatiales (CNES). The reference time period of CNES-CLS18MDT is from 1993 to 2012, and the spatial resolution is 0.125°. Compared with the previous generation of model, this model combines more ocean data through the MOA method [34], which improves the accuracy and resolution of the result, especially in polar and ocean current areas. The reference time periods of these four ocean models are clearly inconsistent, so we first adjusted them to the same time periods using sea level anomaly (SLA) data through the method suggested by Bingham and Haines [35]. Previous studies have successfully used these independent synthetic/ocean datasets to evaluate MDT [21,36,37,38].

3.4. Drifting Buoy Data

Geostrophic velocity was also applied to assess the MDT. The geostrophic velocity was extracted from in situ buoy data. The buoy dataset was provided by the Atlantic Ocean and Meteorological Laboratory (AOML, https://www.aoml.noaa.gov/phod/gdp/index.php (accessed on 28 April 2022)), which was processed by the Kriging method to ensure the quality of the original observations. Buoy data from 1993 to 2012 were used in this study. To obtain the geostrophic velocities from buoy data, the non-geostrophic component must be deducted. Non-geostrophic components in buoy data include Ekman, tidal, inertial, and high-frequency non-geostrophic currents. The Ekman component can be modeled from wind speed and wind stress data. The tidal, inertial, and high-frequency non-geostrophic components can be reduced by a 3-day low-pass filter. After deducting the non-geostrophic components from buoy data, the mean zonal and meridional geostrophic velocities were obtained by averaging the residuals into 0.25° grids. The detained information about obtaining surface geostrophic velocity from buoy data were introduced by Rio et al. [16] and Lumpkin and Johnson [39].

4. Results

4.1. MDT Modeling from the MOA and LS Method

We investigated the performances of the MOA method and the LS method on modeling MDT. Four areas that contain currents were selected as study areas. The DTU21MSS and DIRR6 geoid model were combined to model local MDTs. In the LS method, the observation equation and weight matrix were constructed according to the error information of the GGM and the SNR of the GGM. The geoid was separated into three parts and dealt with accordingly. These three parts represented the geoid signals of three bands, respectively. The first part represented SHs from d/o 2 to 250, the second part from d/o 251 to 300, while the third part expanded SHs from 300 to infinity d/o, where no GGM information was available. As for the construction of the weight matrix, the variance information of the GGM was considered as error models in the LS system, and the smoothness information should be introduced through pseudo-observation considering the SNR of GGM. Moreover, the difference between XGM 2019e_2159 [40] and DIRR6 was used as the diagonal variance information of MSS. Some smoothness information was applied in the observation equation for obtaining a slightly smoother MDT, such as the constraint of MDT gradient norm minimum. It is important to note that the observation equations could be ill-conditioned, due to the large number of SH coefficients. The Schur decomposition [41] was introduced to solve the ill-condition problem. In order to obtain reliable solutions, regularization was carried out, in which the regularization parameters were calculated by L-curve method and the regularization matrix was set as the identity matrix.

The detailed information for modeling MDT by using MOA were shown by Wu et al. [7]. The raw MDT (DTU21 MSS minus DIRR6 geoid) was filtered by a Gaussian filter with a 400 km filter radius to obtain the first guess or initial model. The variance of the MDT was estimated from the residuals (raw MDT minus the CLS18MDT). The variance of the grid point was the variance of the surrounding residuals within a 20° box. Moreover, the a priori error of the MDT was estimated by Bingham’s method. The error was estimated based on the available reference MSS models, GGM models and MDT models. A set of root mean square (RMS) differences between the models we used, and the reference models, were computed as informal error. Then some informal geoid errors and MSS errors were calculated by different reference models, and the error of the MDT was estimated according to error propagation theory. In addition, the second error of the MDT was calculated in the same way as the MSS error. The errors obtained from the reference models of different combinations were compared, and the optimal combination was defined as the smallest RMS of the difference between the two errors. Then the a priori error of MDT was obtained. The covariance of MDT was estimated by CLS18MDT. The correlation radius was the key parameter for estimating the covariance of MDT, which could be determined by fit with an empirical covariance.

4.2. Assessment of MDTs Computed from MOA and LS Method

The MDTs computed by using the MOA (LS) method can be seen in the top (below) panel of Figure 2. The currents can clearly be observed in Figure 2. The patterns of the MDT estimated by the MOA method and the LS method had similar structures. There were some regions that the signals of the MOA-derived MDT were larger than those of the LS-derived MDT, such as the Gulf Current (66°W, 35°N) and Agulhas Current (30°E, 36°S). In order to compare the MDTs estimated by different methods, the synthetic MDT (called ocean data in the following discussion) derived by averaging the ocean reanalysis data and geodetic MDTs models, i.e., SODA, ORAS5, Copernicus, and CNES-CLS18MDT, were used for comparison. Figure 3 shows the discrepancies between the MDTs estimated by the different methods and ocean data.

The differences between the LS-derived MDTs and ocean data were smaller than those of the MOA-derived MDTs. For coastal areas, in the Kuroshio Current area, the differences between the LS-derived MDT and ocean data were clearly smaller than those of the MOA-derived MDT, by a magnitude of about 3 cm, especially for the northern coast of Japan (142°E, 38°N) and southern coast of the Philippines (126°E, 5°N). Notably, for the southern coast of Japan (136°E, 33°N), the discrepancies between the LS-derived MDT and ocean data were larger than those of the MOA-derived MDT. The reason was that this region is located in the Kuroshio Current, where the ocean state is complex, and the ocean data were smooth, which may not be realistic in this area. In the Gulf Current area, the differences between the LS-derived MDT and ocean data were smaller than those of the MOA-derived MDT, by a magnitude of about 4 cm, especially for coast of Cuba (76°W, 17°N). In the Agulhas Current area, compared with the MOA-derived MDT, the main improvement of the LS-derived MDT was located in the eastern and western coasts of Africa, by a magnitude of about 4 cm. In the Greenland Current area, the LS-derived MDT showed smaller differences than the MOA-derived MDT for the southern coast of Greenland (46°W, 60°N), by a magnitude of about 2 cm.

As for ocean current areas, the differences between the LS-derived MDTs and ocean data were slightly larger than for the MOA-derived MDTs. The reason was that the ocean state is complex in ocean current areas, so the ocean data may not be realistic.

Table 1. Statistics of the discrepancies between MDT estimated by different methods and ocean data (units: cm; MOA: multivariate objective analysis; LS: least squares method).

Area	Method	Min	Max	RMS
Kuroshio Current area	MOA	−54.6	15.3	6.2
	LS	−42.4	33.6	4.8
Gulf Current area	MOA	−75.4	40.4	6.9
	LS	−36.0	39.9	4.7
Agulhas Current area	MOA	−59.7	59.6	7.5
	LS	−36.5	32.2	4.9
Greenland Current area	MOA	−58.8	65.3	6.7
	LS	−29.8	28.9	7.2

shows the differences between the MDTs modeled by different methods, and the ocean data. The RMSs of the differences based on the LS method was 4.8 cm, 4.7 cm, 4.9 cm and 7.2 cm for the Kuroshio Current, Gulf Current, Agulhas Current and Greenland Current areas, respectively, which was lower by 1.4 cm, 2.2 cm, 2.6 cm and −0.5 cm, respectively, than for those derived from the MOA-derived MDTs.

Moreover, the values of MDTs along three profiles were extracted for studying the characteristics of MDTs estimated by different methods. Figure 4 shows the profile values that represent the discrepancies between the estimated MDT and the ocean data. The blue lines show the MDT values of the MOA method, and the red lines show the values of the LS method. The profiles are located on the coastal area (profile 1), open sea area (profile 2), and ocean current area (profile 3), which are shown in terms of blue lines in Figure 1. In the coastal area, the profile result shows similar features to the results in Figure 3, where some oscillations appear in these areas, see Figure 4 a(i), b(i), c(i) and d(i). The difference between the LS-derived MDT and the ocean data was smaller than the MOA method. For example, in the Gulf Current area, the results of the LS method were within 5 cm, but the results of the MOA method exceeded 20 cm at 77° W. The results indicated that the MDT derived from the LS method performed better than the MOA method, in the coastal area. In the open sea area, smaller oscillations occurred in the LS-derived results in the Kuroshio Current area and Gulf Current area, see Figure 4 a(ii), b(ii) and c(ii). Spike-like results appeared in the LS-derived profiles much more than in the MOA-derived profiles. The reason was that the ocean data were smooth, and the LS-derived results may preserve more detailed signals. In the ocean current area, the patterns of the LS-derived results and the MOA-derived result were similar. However, the LS-derived results showed larger oscillations and spike-like results. The reason was that in the ocean current area, the sea level change is larger, the sea water flows faster, and the amount of ocean observation data were fewer than for open sea areas, which may lead to a lack of ocean reanalysis data resolution and accuracy; and the ocean data were relatively smooth. It was difficult to distinguish the better results in ocean current area under the complex ocean state in these areas. In terms of the statistics in Table 1, the overall accuracy of the LS-derived MDT was better than the MOA-based MDT.

The better overall accuracy of the LS-derived MDT was because the LS method constructs the design matrix segmentally, based on the error characteristics of the GGM, and then the signals are processed and constrained in different frequency bands to suppress high-frequency noise, which improves the quality of the estimated MDT. The MOA method uses the full available scale signals of geoid (in fact, the signal quality of the GGM is not the same in all frequency bands) and MSS, and the omission errors of the geoid are not handled properly. The error of the geoid is obtained by comparison with four high-degree GGMs. When the input data of MOA method is only the MDT, the MOA method can be seen as an optimal interpolation method. Therefore, the MOA method still has the problem of signal leakage and distortion. Moreover, The LS method is computationally expensive, and takes a long time to calculate. The MOA method can combine other data related to MDT to improve the MDT we estimated.

4.3. Formal Errors of the MDTs Estimated by MOA and LS Method

In the LS method and the MOA method, the MDT error can be calculated by Equation (8) and Equation (11), which can be seen in Figure 5. The formal errors of the LS-derived MDT (in the lower panel of Figure 5) were larger than those estimated by MOA method (in the top panel of Figure 5). The formal errors of the MOA-derived MDT ranged from a few millimeters to 10 cm, with an RMS value of about 3 cm. The formal errors of the LS-derived MDT ranged from about 1 cm to 20 cm, with an RMS value of about 3 cm. In the Kuroshio Current area, the errors of the MDT estimated by the LS method were large in the coastal area, by a magnitude of approximately 20 cm; while the formal errors of the MDT estimated by MOA method were less than 10 cm. The formal errors of the MOA-derived MDT for the southern coast of Japan were larger than for the eastern and southern coastal areas of China, by a magnitude of about 5 cm. The reason was that this region is located in the Kuroshio Current, and the formal errors of the MOA-derived MDT were mainly affected by the current. In the Gulf Current area, the formal errors of the LS-derived MDT were about 15 cm larger than errors of the MOA-derived MDT for the coastal area of North America, and about 20 cm larger than the errors of the coast of Canada. In the Agulhas Current area, the formal errors of the LS-derived MDT were about 15 cm larger than those of the MOA-derived MDT for the coastal area of Africa. In the Greenland Current area, the formal errors of the LS-derived MDT were larger than those of the MOA-derived MDT for the coast of Greenland, by a magnitude of about 20 cm. Moreover, the formal errors of MOA-derived MDT showed large values in ocean current area. Since the error of MSS reached decimeter level over coastal regions, the formal error of MDT computed from the MSS and geoid had at least the same magnitude of MSS’s error through error propagation. Thus, the formal errors of MDT estimated by the LS method may be more reasonable that those derived from the MOA method.

The formal errors of LS-derived MDT and the discrepancies between the estimated MDTs and ocean data (Figure 3) showed high consistency. Both showed larger differences over coastal areas. However, over ocean current areas, different structures could be found between them. The formal errors were relatively small where the discrepancies between the LS-derived MDTs and ocean data were obviously larger, see the structures in (72°W, 38°N) and (18°E, 38°S). The reason was that the error information of MSS we used may not have been accurate enough. In ocean current area, similar forms could be found between the formal errors of MOA-derived MDTs and the discrepancies between the MOA-derived MDTs and ocean data, but different patterns were displayed over coastal areas. The reason was that the variance and a priori error of MDT we used were not accurate, which were smaller than the values in other areas. The results showed that the formal error of the LS-derived MDT may have been more reasonable than that of the MOA-derived MDT.

4.4. Comparison of Geostrophic Velocities Estimated by MDTs Derived from Different Methods

To further evaluate the MDT obtained by different methods, in situ buoy data were introduced as a reference in the form of geostrophic velocity [42]. The geostrophic velocities derived by estimated MDT were filtered by a Gaussian filter with a 60 km filter radius to make them smooth. The differences between the geostrophic velocities calculated by the MOA-derived (LS-derived) MDT and the buoy data are shown in Figure 6 and Figure 7. The LS-derived MDT had better fit with ocean data than the MOA-derived MDT. As for coastal areas, in the Kuroshio Current area, the discrepancies of geostrophic velocities between the LS-derived results and buoy data were clearly smaller than the MOA-derived results. For zonal geostrophic velocities, the discrepancies of geostrophic velocities between the LS-derived results and buoy data were about 2 cm/s smaller than that of the MOA-derived results for the coast of the Philippines. For meridian geostrophic velocities, the discrepancies of the LS-derived results were 2 cm/s smaller than the MOA-derived results for the coast of Japan. In the Gulf Current area, the discrepancies of geostrophic velocities between the LS-derived results and buoy data were clearly smaller than the MOA-derived results for the coasts of Canada and Cuba, by a magnitude of 3 cm/s. In the Agulhas Current area, the discrepancies of geostrophic velocities between the LS-derived results and buoy data were smaller than the MOA-derived results for the eastern and western coasts of Africa, by a magnitude of 3 cm/s. In the Greenland Current area, the discrepancies of geostrophic velocities between the LS-derived results and buoy data showed few discrepancies from that of the MOA-derived results. The reason was that there are several currents in the area, the sea states are complicated, and there were less buoy data than other areas, which may have resulted in inaccurate geostrophic velocity derived from buoy data. In addition, we have mentioned that the difference between the LS-derived MDT and the MOA-derived MDT was small.

As for ocean current areas, the discrepancies of geostrophic velocities between the LS-derived results and ocean data were slightly larger than the MOA-derived results, such as in the Gulf Current. The reason may have been because the sea surface state changes rapidly in the ocean current region, and the accuracy and resolution of buoy data were not high enough. As for open sea areas, the differences of geostrophic velocities between the LS-derived results and ocean data were similar to the MOA-derived results, which ranged from about −1.5 cm/s to 1.5 cm/s in the Kuroshio Current area, from −1 cm/s to 1 cm/s in the Gulf Current area, from −2 cm/s to 2 cm/s in the Agulhas Current area, and from −1 cm/s to 1.5 cm/s in the Greenland Current area.

The statistics of the discrepancies between the geostrophic velocities computed by the MDTs based on the different methods and buoy data are shown in

Table 2. Statistics of the discrepancies between the geostrophic velocities computed from different MDTs based on different methods and buoy data (units: cm/s; u: zonal velocities; v: meridian velocities; MOA: multivariate objective analysis; LS: least squares method).

Study Area	Method	Geostrophic Velocities	Min	Max	RMS
Kuroshio Current area	MOA	u	−104.7	99.0	16.1
		v	−103.7	169.1	15.4
	LS	u	−116.6	75.1	15.7
		v	−98.9	95.0	13.4
Gulf Current area	MOA	u	−66.9	51.9	10.1
		v	−156.1	51.7	12.1
	LS	u	−63.6	50.4	9.9
		v	−151.2	61.1	11.8
Agulhas Current area	MOA	u	−126.3	110.9	19.7
		v	−171.8	129.2	23.3
	LS	u	−124.9	113.3	18.6
		v	−119.0	123.6	21.0
Greenland Current area	MOA	u	−77.5	48.6	8.4
		v	−58.7	68.6	9.0
	LS	u	−77.5	48.6	8.4
		v	−58.7	68.6	9.0

. The RMSs of the discrepancies between the zonal (meridional) velocities calculated by LS-derived MDT and the buoy data were 0.4 cm/s (2.0 cm/s), 0.2 cm/s (0.3 cm/s) and 1.1 cm/s (2.3 cm/s) over the Kuroshio Current area, the Gulf Current area and the Agulhas Current area, respectively, which were lower than the velocities calculated by MOA-derived MDT. Meanwhile, the LS-derived MDT and the MOA-derived MDT had comparable performances for the Greenland Current area, as the RMSs of the discrepancies between the geostrophic velocities calculated by these two MDTs and the buoy data were 8.4 cm/s in zonal and 9.0 cm/s in meridional. These results indicated that the geostrophic velocity derived from LS-derived MDT outperformed that of MOA-derived MDT, especially in coastal and ocean current areas.

5. Discussion

The other MDT model and ocean models were applied to evaluate the MDT obtained by different methods. CNES-CLS18MDT was the reference model, which was modeled by combining the MSS model, geoid model, buoy data and hydrological profiles data. The differences between the geostrophic velocities calculated by the MOA-derived (LS-derived) MDT and those of CNES-CLS18MDT are shown in Figure 8 and Figure 9. The LS-derived MDT had better fit with the CNES-CLS18MDT than the MOA-derived MDT. As for coastal areas, in the Kuroshio Current area, the discrepancies of geostrophic velocities between the LS-derived results and buoy data were clearly smaller than those of the MOA-derived results. The discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were about 10 cm/s smaller than those of the MOA-derived results for the coasts of the Philippines and Japan. In the Gulf Current area, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were similar. In the Agulhas Current area, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were smaller than those of the MOA-derived results in the eastern and western coasts of Africa, by a magnitude of 20 cm/s. In the Greenland Current area, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results showed few discrepancies from those of the MOA-derived results. For ocean current areas, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were slightly larger than those of the MOA-derived results, such as in the Gulf Current. For open sea areas, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were similar to those of the MOA-derived results. These results were similar to the results mentioned in Section 4.4.

The statistics of the discrepancies between the geostrophic velocities computed by the MDTs based on different methods and those of the CNES-CLS18MDT are shown in

Table 3. Statistics of the discrepancies between the geostrophic velocities computed from different MDTs based on different methods and that of CNES-CLS18MDT (units: cm/s; u: zonal velocities; v: meridian velocities; MOA: multivariate objective analysis; LS: least squares method).

Study Area	Method	Geostrophic Velocities	Min	Max	RMS
Kuroshio Current area	MOA	u	−95.8	123.7	9.7
		v	−127.1	115.4	10.9
	LS	u	−45.3	74.7	5.9
		v	−54.6	48.9	5.9
Gulf Current area	MOA	u	−29.3	26.2	4.3
		v	−72.1	22.3	4.5
	LS	u	−21.4	13.9	2.7
		v	−68.9	16.3	3.8
Agulhas Current area	MOA	u	−79.1	140.2	10.7
		v	−176.1	125.3	14.0
	LS	u	−14.9	36.3	4.6
		v	−28.2	29.3	3.7
Greenland Current area	MOA	u	−26.1	14.6	2.9
		v	−29.6	26.4	2.7
	LS	u	−19.7	12.6	2.7
		v	−12.2	15.8	2.4

. The RMSs of the discrepancies between the zonal (meridional) velocities calculated by LS-derived MDT and the CNES-CLS18MDT were 3.8 cm/s (5.0 cm/s), 1.6 cm/s (0.7 cm/s), 6.1 cm/s (10.3 cm/s) and 0.2 cm/s (0.3 cm/s), over the Kuroshio Current area, the Gulf Current area, the Agulhas Current area and the Greenland Current area, respectively, which were lower than that of the velocities calculated by MOA-derived MDT. These results indicate that the geostrophic velocity derived from LS-derived MDT outperforms that of MOA-derived MDT, especially in coastal area and ocean current areas.

Moreover, the ocean reanalysis models, such as SODA, ORAS5 and Copernicus, which collect the geostrophic velocities data, were applied to evaluate the MDT obtained by different methods. In order to reduce the systematic error between the ocean reanalysis models, the mean of the three ocean reanalysis models was set as the reference data (called Ref-ocean-model in the following discussion). The differences between the geostrophic velocities calculated by the MOA-derived (LS-derived) MDT and the Ref-ocean-model are shown in Figure 10 and Figure 11. The geostrophic velocities derived by LS-derived MDT had better fit with the Ref-ocean-model than the MOA-derived MDT. As for coastal areas, the discrepancies of geostrophic velocities between the LS-derived results and Ref-ocean-model were mostly smaller than those of the MOA-derived results, especially for the coasts of the Philippines, Japan, and eastern and western coasts of Africa. For ocean current areas, the discrepancies of geostrophic velocities between the LS-derived results and the Ref-ocean-model were slightly larger than those of the MOA-derived results, such as in the Gulf Current. For open sea areas, the discrepancies of geostrophic velocities between the LS-derived results and the CNES-CLS18MDT-derived results were similar to those of the MOA-derived results. These results were similar with the results mentioned above, which indicated that the conclusion was not affected by the reference model.

The statistics of the discrepancies between the geostrophic velocities computed by the MDTs based on different methods and the Ref-ocean-model are shown in

Table 4. Statistics of the discrepancies between the geostrophic velocities computed from different MDTs based on different methods and the Ref-ocean-model (units: cm/s; u: zonal velocities; v: meridian velocities; MOA: multivariate objective analysis; LS: least squares method).

Study Area	Method	Geostrophic Velocities	Min	Max	RMS
Kuroshio Current area	MOA	u	−106.4	151.5	12.2
		v	−133.7	155.1	13.4
	LS	u	−59.0	74.9	9.1
		v	−64.2	99.3	8.9
Gulf Current area	MOA	u	−60.7	33.5	5.4
		v	−108.0	33.0	7.2
	LS	u	−67.1	33.6	5.3
		v	−96.7	41.2	6.7
Agulhas Current area	MOA	u	−77.2	136.2	12.4
		v	−192.7	121.7	15.0
	LS	u	−29.2	59.6	7.5
		v	−58.3	52.9	6.9
Greenland Current area	MOA	u	−28.0	23.8	3.9
		v	−30.3	26.7	4.1
	LS	u	−28.1	26.7	3.7
		v	−19.3	22.5	3.9

. The RMSs of the discrepancies between the zonal (meridional) velocities calculated by LS-derived MDT and the CNES-CLS18MDT were 3.1 cm/s (4.5 cm/s), 0.1 cm/s (0.5 cm/s), 4.9 cm/s (8.1 cm/s) and 0.2 cm/s (0.2 cm/s) over the Kuroshio Current area, Gulf Current area, Agulhas Current area and Greenland Current area, respectively, which were lower than the velocities calculated by MOA-derived MDT. These results indicated that the geostrophic velocity derived from LS-derived MDT outperformed that of MOA-derived MDT, especially in coastal and ocean current areas.

6. Conclusions

We focused on the comparison of methods of modeling MDT by using the multivariate objective analysis (MOA) method and rigorous least squares (LS) method. In particular, we evaluated the applicability of these two methods and compared their performances over different oceanic areas with different ocean state and hydrological conditions. Four local MDTs were computed, and the estimated MDTs were assessed by independent ocean data and buoy data. Moreover, the formal errors of the estimated MDTs based on these two methods were also analyzed and compared. The numerical results showed that:

(1)

The MDT derived from the LS method outperformed the MDT computed from the MOA method, especially over coastal areas and ocean current areas. The RMSs of the discrepancies between the LS-derived MDT and ocean data were 4.8 cm, 4.7 cm, 4.9 cm and 7.2 cm, for the Kuroshio Current area, Gulf Current area, Agulhas Current area and Greenland Current area, respectively, which were lower than those of the MOA-derived MDT, by a magnitude of 1.4 cm, 2.2 cm, 2.6 cm and −0.5 cm, respectively. The reason is that the LS method constructs the design matrix segmentally based on the error characteristics of the GGM, and then the signals are processed and constrained in different frequency bands to suppress high-frequency noise, which improves the quality of the estimated MDT;

(2)

The formal error of the MDT estimated by the LS method was more reasonable than that estimated by the MOA method. The errors of the MDT estimated by the LS method were prominent over coastal areas, which have larger magnitude, than estimated by the MOA method. The patterns of the formal errors of the LS-derived MDT were more realistic, since the errors of MSS models usually exceeded decimeter level along the coast, indicating the formal error of MDT computed from the MSS and geoid has at least the same magnitude of error as MSS, through error propagation;

(3)

Moreover, the geostrophic velocity derived from the LS-derived MDT was better than from the MOA-derived MDT, especially over coastal regions and ocean current areas. The RMSs of the discrepancies between the zonal (meridional) velocities calculated by the LS-derived MDT and the buoy data were 0.4 cm/s (2.0 cm/s), 0.2 cm/s (0.3 cm/s) and 1.1 cm/s (2.3 cm/s) smaller than of the velocities calculated by the MOA-derived MDT over the Kuroshio Current area, Gulf Current area and Agulhas Current area, respectively. The comparison between geostrophic velocities estimated by MDTs derived from different methods and the ocean models showed similar results. The results indicate that the LS-derived MDT outperforms the MOA-derived MDT.