Spectrally Consistent Mean Dynamic Topography by Combining Mean Sea Surface and Global Geopotential Model Through a Least Squares-Based Approach

The filtering procedure is usually mandatory for modeling mean dynamic topography (MDT) when a geodetic approach based on the Mean Sea Surface (MSS) and the Global Geopotential Model (GGM) is used. This is due to the inconsistent spectral contents between MSS and GGM. However, traditional isotropic filtering algorithms (e.g., Gaussian filter) consider neither the MDT locations nor their azimuth when smoothing the signal within the filtering radius. Hence, the isotropic filtering will attenuate the MDT signal near the current and filter the current signal into the surrounding ocean, which may lead to signal contamination and distortion. In this study, we set up a least squares-based (LS) approach to model MDT signal from the altimeter-derived MSS and geoid height using spherical harmonics from GGMs, where MDT is parameterized by Lagrange Basis Functions (LBFs). The design matrix is segmentally established, considering the error information of GGM in various spectral bands. Numerical experiments in the Gulf Stream show that applications of full error variance-covariance matrix or only diagonal error variance of GGM may have marginal effects on the MDT modeling. The MDT computed from this LS-based approach using the latest releases of Gravity Field and Steady-state Ocean Circulation Explorer (GOCE) geoid models, i.e., GO_CONS_GCF_2_DIR_R6 and Gravity Observation Combination 06s model (GOCO06s), have the best agreement with the comparison data, especially near the current region. Deduced geostrophic velocities based on the MDT solutions show that the LS-based approach recovers the current signal better than the Gaussian filtering by 1.8 cm/s. Estimated error map illustrates that errors are more concentrated near the coastal region.


INTRODUCTION
The mean dynamic topography (MDT) plays a very important role in ocean circulation, global climate change, and vertical datum unification (Le Traon et al., 2015;Woodworth et al., 2015). In the geodetic sense, the MDT is represented as the difference between the mean sea surface (MSS) and the geoid in a specific reference period, which emphasizes the importance of the accurate knowledge of MSS and geoid (Rio, 2004;Andersen and Knudsen, 2009;Knudsen et al., 2011).
Since the first operation of the altimetry satellite in 1970s, there have been more than 20 altimetry missions used for global sea level monitoring, giving a centimeter-level accuracy in the open ocean (Vignudelli et al., 2019;Donlon et al., 2021). The accumulation of satellite altimetry data over 20 years greatly refined the MSS (Hwang et al., 2002;Andersen and Knudsen, 2009). Consequently, the latest models such as Centre National d'Etudes Spatiales -Collect Localisation Satellites 2015 Mean Sea Surface model (CNES-CLS15MSS) (Schaeffer et al., 2016) and Danmarks Tekniske Universitet 2018 Mean Sea Surface model (DTU18MSS) (Andersen et al., 2018b), combining 20 years of altimetric data, can refine MSS at the spatial resolution down to 10 km. With the operation of the Gravity Recovery And Climate Experiment (GRACE) mission (Swenson and Wahr, 2002;Tapley et al., 2004) launched in 2002 and GRACE Follow-on (GRACE-FO) mission launched in 2018 Landerer et al., 2020), especially the operation of the subsequent Gravity field and steady-state Ocean Circulation Explorer (GOCE) launched in 2009 (Drinkwater et al., 2003;Bock et al., 2014;Brockmann, 2014;Wu et al., 2017), the modeling accuracy and spatial resolution of the global gravity field have been improved unprecedentedly (Pail et al., 2010;Bruinsma et al., 2014;Tziavos et al., 2015). The GRACE and GRACE-FO missions have an initial orbit altitude of about 500 km and can accurately measure the long-wavelength of global gravity field signal, which made it possible for the first time to obtain reasonable MDT results relying solely on satellite gravity data (Vianna et al., 2007;Knudsen et al., 2011). In 2009, the European Space Agency (ESA) launched a GOCE mission that carried a gravity gradient payload for the first time (Drinkwater et al., 2003;Baur and Grafarend, 2006). With an orbit altitude of onlỹ 250 km, the GOCE mission can sense higher-frequency gravity field signal up to degree and order (d/o)~250. Currently, the maximum expansion of the Global Geopotential Model (GGM) based on pure GRACE gravity data is d/o~200 (Tapley et al., 2003;Mayer-Gürr et al., 2014;Kvas et al., 2019a), and due to the complementarity of gravity signal from the two gravity satellite missions, the pure satellite GGM integrating data from GRACE/GOCE can be expanded to d/o 300, corresponding to a spatial resolution of about 70 km.
The spatial resolution of GGM is still much lower than that of MSS (Bingham et al., 2014). MSS is furthermore given as a uniform grid over ocean, while GGM needs to be represented as the sum of global Spherical Harmonics (SHs), so the spectral content between MSS and GGM will not be completely matched (Bingham et al., 2008). Deriving the MDT by subtracting the two models will inevitably cause signal leakage problems. To resolve this, several researchers have utilized filtering in either space or spectral domains to obtain smooth MDT solutions by removing high-frequency information. The Gaussian filter, which is essentially an isotropic weighted low-pass filter, was first proposed by Jekeli (1981). It has been widely used in MDT modeling and the choice of filtering parameters has been discussed in detail (Knudsen et al., 2011;Siegismund, 2020). To avoid the signal contamination and attenuation caused by the isotropic filters, Bingham (2010) and Sánchez-Reales (2016) represent the filter as an anisotropic diffusion process reducing the MDT attenuation near currents.
In order to accurately determine a geodetic MDT and comprehensively analyze its error level, it is necessary to take into account the error information of input data in the MDT modeling. On the one hand, the commission error increases significantly with the expansion of GGM, which should be assessed and reweighted in the MDT modeling based on the error information of GGM. On the other hand, the error models of MSS and GGM could help to evaluate the formal errors of the MDT. However, the evaluation of the error characteristic of the MDT would still be hampered due to the filtering procedure; hence, methods that could evaluate errors of the MDT are presented in some studies (Bingham et al., 2014;Bai et al., 2020). For instance, Rio et al. (2011)and Rio et al. (2014) introduced an objective analysis method considering the variance-covariance and the correlation between observations for MDT modeling, which helps suppress the noise level in the MDT solution. Alternatively, Becker et al. (2012)and Becker et al. (2014) established a complete observation equation considering the error characteristics of MSS and GGM. As a result, a smooth MDT solution and its associated errors were estimated without any additional filtering procedure (Chambers et al., 2017).
Inspired by this approach, this study set up a Least Squares-based (LS) approach to combine MSS and GGM for MDT estimation, which avoids additional smoothing filters after the LS adjustment. The satellite-only GGMs are used for MDT modeling. Both the omission error and the commission error information of the GGM are considered segmentally and introduced into observation equations, so that a spectrally consistent MDT could be estimated in the LS system. The MDT solutions using full error variance-covariance information or diagonal error variance information of GGM in the MDT modeling are also compared. The effects of different choices of GGMs are investigated quantitively, and finally geostrophic velocities are validated. This paper starts with a brief introduction of the principle of the LS-based approach for MDT modeling. Section 2 introduces the research area and data used in this study. The construction of the complete observation equations and the weight matrices based on error information provided by input models is detailed in Section 3. Section 4 illustrates estimated MDT solutions based on the LS-based approach, where the effects of different choices of GGMs are analyzed. Deduced geostrophic current velocities are validated, subsequently. Finally, a brief conclusion is drawn in Section 5.

Research Area
Our research area is the Gulf Stream (15°N-55°N, 40°W-80°W) illustrated in Figure 1 [image is adapted from the General Bathymetric Chart of Oceans (GEBCO) world map 2014, www.gebco.net], with the distribution of ocean currents. The Gulf Stream is one of the strongest ocean currents in the world and has a significant influence on climate, mass/heat transportation, and ocean fishery (Rossby et al., 2010;Palter, 2015). Formed by the Coriolis force and the trade wind, the Gulf Stream originates from the North equatorial current and is obstructed by the American continent, which flows northward first and then eastward through the Gulf of Mexico, where the maximum speed could reach 2 m/s. When flowing in deep water, a large number of currents shed off from the main current, forming mesoscale eddies (Kang and Curchitser, 2013). As a result, the variation of the MDT around the Gulf Stream is complex, leading to a difficulty in the MDT signal recovery from the geodetic approach (Rio, 2004;Klymak et al., 2016). On the other hand, the complexity of the currents in this area provides favorable conditions for validating the MDT modeling algorithm, especially its performance near current and coast, leading to the Gulf Stream becoming a hot spot in MDT research.

MSS and GGM
In order to accurately separate the MDT signal from MSS and GGM, the latest MSS of CNES-CLS15MSS is introduced for modeling (Schaeffer et al., 2016;Pujol et al., 2018), where 20 years of altimetric data from 1993 to 2012 are assimilated. The CNES-CLS15MSS is represented on a 1/8°spatial resolution grid. The reference time period is 1993-2012, and the reference ellipsoid is the World Geodetic System 1984 (WGS84).
The ITSG-Grace2014s GGM (ITSG), published by the Graz University of Technology, provides full error variance-covariance information that can be conducive to the investigation of the commission error modeling (Mayer-Gürr et al., 2014). However, the ITSG is a GRACE-only GGM with a maximum expansion of d/o 200, which will be a limitation in the MDT modeling, spectrally ( Table 1). In the following the ITSG model will be called the "GRACE-only GGM." In addition, with the completion of a GOCE mission, the accuracy of GGM based on satellite-only gravimetric data has been unprecedently improved, especially in the open ocean. Therefore, the Gravity Observation Combination 06s model (GOCO06s) that provides full error variancecovariance matrix is also introduced for MDT modeling (Kvas et al., 2019b).
ESA has released six releases of direct GGM solutions based on GRACE/GOCE gravity data, where the fourth-release (GO_CONS_GCF_2_ DIR_R4, DIR4), the fifth-release (GO_CONS_GCF_2_DIR_R5, DIR5), and the sixth-release (GO_CONS_GCF_2_DIR_R6, DIR6) are introduced for MDT modeling. Notably, the DIR6 recalibrates the satellite orbits compared to the DIR5 (Förste et al., 2019). Besides, in order to better assess the performance of the LS-based approach in the MDT estimation, the timewise solutions that only assimilated GOCE data, including the fifth-release (GO_CONS_GCF_2_TIM_R5, TIM5) and the sixth-release (GO_CONS_GCF_2_TIM_R6, TIM6, and GO_CONS_GCF_ 2_TIM_R6e, TIM6e), are also introduced for MDT estimation. In the following the releases of GOCE GGM will be called the "GOCE-based GGM." Table 1 summarizes more detailed information about these GGMs. In addition, the XGM 2019e_2159 combined GGM (Zingerle et al., 2020), provided by Technical University of Munich, is introduced as the complementary information in the weight matrix construction in the LS-based approach. The reason for choosing XGM 2019e_2159 is that the XGM 2019e_2159 assimilated the latest Altimetry data, satellite gravity data, and ground data for modeling, which could be expanded up to 2190d/ o. And comparing with the other high d/o GGMs e.g., EIGEN6c4  or EGM 2008 (Pavlis et al., 2012), the XGM 2019e_2159 shows a better performance in the ocean, especially in the coastal regions (Zingerle et al., 2020).

Synthetic MDT
In order to assess the performance of estimated MDT, the latest Centre National d'Etudes Spatiales -Collect Localisation Satellites 2018 Mean Dynamic Topography model (CNES-CLS18MDT, CLS18) and the Danmarks Tekniske Universitet 2017 Mean Dynamic Topography model (DTUc17MDT, DTU17) with their associated geostrophic velocity fields are introduced for comparison and validation (Knudsen et al., 2021;Mulet et al., 2021). The CLS18 utilized CNES-CLS15MSS and GOCO05S for modeling. All drifting buoy velocities and hydrological profiles from 1993 to 2017 are assimilated in CLS18 for current signal augmentation. Estimated error of the CNES-CLS18 could reach 10 cm around current areas and down to~2 cm in the other areas (Mulet et al., 2021). The DTU17 has been developed using the DTU15MSS as MSS, and EIGEN-6C4 for geoid modeling. Drifter data from 1992 to 2015 are assimilated in the DTU17, and the quasi-gaussian filter (Knudsen et al., 2021) was implemented to best fit the velocities of oceanographic drifting buoys. The estimated errors of the DTU17cMDT range from a few centimeters to a few decimeters, and they exceed 20 cm in areas of the major ocean currents (Knudsen et al., 2021). The reference period of both CLS18 and DTU17 is 1993-2012, with the spatial resolution of 1/8°. Table 2 shows more details about these two models. In this study, we take geostrophic velocity fields provided by CLS18 and DTU17 as validation data, considering that both models assimilated a large amount of in-situ data. A synthetic MDT model (Synthetic_MDT) and associated geostrophic velocity field (Synthetic_Current) were calculated from the average of the CLS18 and DTU17 for comparison and validation.

Gaussian Filtered MDT
The Gaussian filtered MDT will also be introduced for comparison in this study, which is derived based on a direct solution combining GGM and CNES-CLS15MSS. In this study, the same GGM information are introduced to estimated MDT and the Gaussian filtered MDT to ensure that the content between the two MDT solutions is comparable. The Gaussian filter used in this study is introduced from Jekeli (1981) and Wahr et al. (1998). The MDT values (X) could be Gaussian filtered bŷ whereX i (i 1, 2, 3, ...) represents the filtered MDT value and W represents the weighting vectors. The weighting factors read based on Wahr et al. (1998) as where α and R filter represents the spherical distance of the MDT point and the filtering radius (m). The filtering radius (R filter ) for the Gaussian filtered MDT is given as (Barthelmes, 2013): where l max denotes the max degree of the expansion of SHs.

METHODOLOGY
We initially set up an LS-based approach to separate MDT signal from MSS and GGM. The key point is the establishment of a complete observation equation and the associated weight matrix, which will prominently influence the separation of MDT signal from MSS and GGM. Here, we express GGM and MDT in the spatial domain, which is given as (Bingham et al., 2008) where θ and λ represent the latitude and longitude in spherical coordinate system, respectively. Equation 5 illustrates the model used for MDT reconstruction, and the geoid and the MDT should be parameterized in a uniform grid. Consequently, the geoid is expanded based on a series of SHs, which is given as where GM represents the gravitational constant times the Earth's mass, R denotes the Earth's radius, γ(B) denotes the normal gravity at geodetic latitude B, P nm (cos θ) denotes the fully normalized Legendre polynomial at degree n and order m, and C nm , S nm is the corresponding Stokes coefficients. In addition, the MDT is parameterized by the Lagrange Basis Function (LBF) (Becker et al., 2012): where K represents the number of basis function b k and a k represents the MDT values at nodes ( θ, λ). Due to the computer limitation, LBF with 4 parameters was used to interpolate the MDT in this study, despite the recommendation to use 16 parameters in Shi et al. (2020). MSS and GGM are extracted on a 0.5°grid, respectively, in order to reduce the correlation between the grids. For each unknown MDT point, the four points around it in the 0.5°grid are found, and a local coordinate system in this 4-points-rectangle is established. Then, the unknown MDT value in this local system is calculated based on the interpolation using LBF with 4 parameters. Finally, all the unknown MDT values can be parameterized by LBF, and the parameter coefficient matrix (A mdt ) is derived. More details of the MDT parameterization can be seen in Supplementary Appendix SA, Becker et al. (2012), and Shi et al. (2020). Equation 5 can be rewritten as where v represents the residuals derived from the LS system, and x cs and x mdt denote the unknown SH coefficients and the MDT values, which should be estimated based on the LS theory. The determination of A cs and A mdt can be seen in Supplementary Appendix SB.
In this study, we investigate the use of the satellite-only GGM. In principle, SHs from d/o 2 to ∞ should be used for expression of the geoid. However, the Signal to Noise Ratio (SNR) of GGM will decrease with the increase of the expansion of GGM (Tsoulis and Patlakis, 2013;Bruinsma et al., 2014). Consequently, we segmentally construct observation equation and the weight matrix based on the error information provided by GGMs.
The expression of GGM is split into three parts and processed accordingly. The first parameters group of SHs (cs1) represents SHs from d/o 2 to an appropriate cut-off frequency of GGM (e.g., where SNR>1), the second group (cs2) serves as a buffer between cs1 and cs3, and the third group (cs3) where no GGM information is available expands SHs from max d/o of cs2 to infinity d/o, hence: In order to better extract the MDT signal from MSS and GGM, more smoothing conditions and a priori information should be added, which will be discussed in the following sub-sections.

The Construction of the Weight Matrix
We directly introduce SH coefficients C nm , S nm of GGM as a priori information for cs1, which is given as Here, the variance information provided in GGM (K cs1 ) is considered as error models in the LS system. Accordingly, the commission error introduced by GGM would be reweighted and suppressed in the LS system.
The smoothness information should be introduced through pseudo-observation considering the SNR of GGM. Supposing that the random variables are normally distributed, the degreewise variance function of the SH coefficients can be expressed based on Kaula's rule (Kaula, 1966): where σ 2 n represents the variance information at degree n. The smoothness information is introduced in cs2 and cs3, and the variance matrix results in Consequently, the pseudo-observation equations can be expressed as In the combination of MSS and GGM, the main obstruction is the limited spatial resolution of GGM. Consequently, information from cs3 is usually neglected or set to zero, and a spatial or spherical filter is introduced for smoothing (Knudsen et al., 2011;Ke et al., 2019). Hence, considering the omission errors of GGM, the combined MDT would suffer from a serious signal leakage problem. As a result, it is important to consider the error information for cs3 for a spectral consistent MDT.
In our approach, cs3 is treated separately (Eq. 9), where S A cs3 · x cs3 is the signal from cs3, thus, where MSS MSS − S. In this study, the discrepancy between XGM 2019e_2159 (Zingerle et al., 2020) and GGM used in the LS-based approach as the diagonal variance information of MSS (K MSS ).

The Additional Smoothness Information
In order to derive a slightly smoother MDT and associated geostrophic current velocity field, additional smoothness information needs to be supplemented in the observation equation. In this study, we try to minimize the norm of the gradient of the MDT (Becker et al., 2014). The parameterization matrix of the MDT part in Eq. 8 can be expressed as where A mdt denotes the parameterization matrix based on LBF. The derivative of the parameterized MDT in zonal (∇A x ) and meridian (∇A y ) direction can be expressed as Hence, the smoothing information that could reduce the norm of the MDT gradient could be supplemented in observation equation by

The Complete Observation Equation
This study expresses MSS as the combination of the MDT and the geoid, where the geoid is represented as a sum of SH functions and the associated error information provided by GGM is introduced in the modeling. The MDT is parameterized and estimated in a LS system, and observation equations are segmentally established considering the SNR of GGM. It is notable that Eq. 9 may suffer from ill-condition problem, since the SH coefficients are introduced in observation equations. Also, the estimated MDT using the LS-based approach may be coarse, due to the noises introduced from the MDT parameterization and the weight matrix. Therefore, the SH coefficients of the satellite-only GGM (Eq. 10) and the gradient term of the MDT (Eq. 19) are introduced in this study as constraints to help derive a reasonable MDT solution. In the end, the complete observation equation with its weight matrix can be expressed as Therefore, observation equation can be solved by the weighted LS method:x wherex represents the unknown parameters in Eq. 20, and A, P represent the design matrix and the weight matrix, respectively. However, this study focuses on the estimation of the MDT parameters (x mdt ), so the calculation of Eq. 22 is time consuming and may lead to ill-condition problems. Hence, the Schur decomposition (Lawson and Hanson, 1995) is introduced to help solve the solution. Supplementary Appendix SC shows more information about the derivation of the x mdt . As a result, this Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 795935 6 approach avoids introduction of additional filtering procedure after the LS adjustment and is expected to retain more detailed signal. Further, the posterior variance of the MDT solution can also be estimated by error propagation. This is beneficial to further assimilation by combining information such as ocean numerical data (Freiwald, 2013;Rio et al., 2014).

RESULTS
The performance of the LS-based MDT modeling approach is evaluated in the Gulf Stream area (Figure 1). We first investigate the importance of the full variance-covariance information of GGM in the MDT modeling, where GRACE-only GGM and GOCE-based GGM are used. Then, we study the effect of gravity signal in the modeling of MDT using the latest seven releases of GOCE-based GGMs, and estimated errors of the MDT solutions will be analyzed and assessed. Effects of different choices of the cut-off frequency for cs1 will be investigated.

Importance of the Full Error Variance-Covariance Matrix
As for GRACE-only GGM (ITSG), the choice of the cut-off frequency for segmentation is cs1: d/o 2-160, cs2: d/o 161-180, while the choice of the cut-off frequency for GOCE- Considering the cut-off frequency of cs1, the radius of the Gaussian filtered MDT is chosen as 1.1°(~125 km) for ITSG, and 0.65°(~75 km) for GOCO06s, based on Eq. 4. The MDT solutions with 1°resolution are estimated and compared in the following experiment.
Full error variance-covariance matrix and diagonal error variance matrix from ITSG and GOCO06s are introduced in the LS system, respectively, and the corresponding solutions are compared. The MDT solutions estimated using ITSG are shown in Figure 2. Figure Table 3. Notice that the land area is ignored in this experiment.
All MDTs provide reasonable solutions, with similar patterns of the signal distribution. Specifically, bounded by the Gulf Stream (Figure 1), the MDT in its southern part is above 0.4 m, while the MDT in the northern part is below 0.4 m. Also, it can be seen from Figure 2A, that the variation of the MDT around (38°N/70°W) is more than 1 m, and we will take it as a comparison point in this study. The MDTs derived by the two approaches could reach the level of -0.8-1.2 m in the research area. The Gaussian filtered MDT is apparently smoother than estimated MDT, while there is sharper signal in the estimated solution. The reason is that the Gaussian filter operates an isotropic spatial weighting average procedure to the signal for a smoother solution (Bingham, 2010). On the contrary, there are no additional filters introduced in estimated MDT after the LS adjustment.
As can be seen from Figures 2E-G, misfits of all three MDT solutions can reach the level of ±0.5 m. However, their distribution is different. Misfits between the Gaussian filtered MDT and the Synthetic_MDT in Figure 2E show a clear signal leakage, especially near the current and the coastal area. Specifically, the signal near current is weakened for about 0.3 m. This is because the isotropic filtering causes attenuation of the MDT gradient signal, resulting in a large discrepancy between the MDT solution and the Synthetic_MDT. The loss of the signal will have an influence not only on the model accuracy but also on further investigation in ocean currents.
In addition, Figures 3A-C show misfits between Gaussian filtered MDT, estimated MDT with full error information, and diagonal error information of GOCO06s, against the Synthetic_MDT, respectively. The statistics are shown in Table 4. Apparently, the improvement in GOCE-based MDTs is significant compared with the GRACE-only MDT, where the SD of misfits decreases from 0.046 m (ITSG) to 0.037 m (GOCO06s). It is seen from Figure 3B that there are highfrequency misfits along 38°N in the estimated MDT, but it's apparently better suppressed comparing with that in Figure 2F using the GRACE-only GGM, which is the contribution of GOCE signal.
However, both results obtained from ITSG and GOCO06s show that there are only marginal differences between the MDT solutions using full error variance-covariance of GGM and using only diagonal error variance information for error modeling. Statistics show a maximum difference of 0.6 cm between these two MDT solutions, with a SD of~0.05 cm. This is likely because the diagonal elements of SHs coefficients play a major role in GGM signal (Gruber, 2001;Ran et al., 2021). Consequently, only the diagonal error variance information of GGM will be considered in the following MDT modeling using GOCEbased GGM in the next section.
Comparing with the Synthetic_MDT, both the Gaussian filtered solution and estimated MDT solution will show signal attenuation near the current. For example, alternating positive and negative misfits are seen around 30°N-38°N, which is induced by the over-smooth in the Gaussian filtering or the smoothness condition in the LS system. The estimated MDT recovers more signal than the Gaussian filtered one. Misfits of the estimated MDTs ( Figures 2F, 3B) around 30-38°N are significantly smaller compared with the Gaussian filtered solution (Figures 2E, 3A). Tables 3, 4 show that, although the SD of misfits between Gaussian filtered MDT and the Synthetic_MDT decreases from 0.084 m (ITSG) to 0.046 m (GOCO06s), the MDT estimated by the LS-based approach are still a better fit with the Synthetic_MDT. The SD of misfits between the estimated MDTs and the Synthetic_MDT decreases from 0.046 m (ITSG) to 0.037 m (GOCO06s). In addition, the signals around the coast in the southwestern part of the research area are better recovered than the Gaussian filtered solution, indicating the superiority of the LS-based approach in the MDT modeling.

Improved GOCE-Based GGMs
In order to further assess the performance of GOCE-based GGMs in MDT estimation, seven GOCE-based GGMs (i.e., DIR4, DIR5/ TIM5, DIR6/TIM6/TIM6e, and GOCO06s) are investigated, and their MDT solutions are compared. The strategies for MDT modeling are the same as in Section 4.1. Considering the SNR of the GGMs, the cut-off frequency for cs1 is chosen as d/o 2-260, and cs2 is chosen as d/o 261-280. Again, the signal using the difference between GOCE-based GGMs and XGM 2019e_2159 are treated as the variance information of MSS in observation    Table 5. The Gaussian filtered MDT still suffers from the severe signal leakage problem near the current and coast region in the southwestern part of the research area. For instance, misfit around 38°N/70°W is~±0.1 m in the estimated solutions, while it increases to~±0.2 m in the Gaussian filtered solution. The statistics in Table 5 show that the SD of misfits of the MDT solutions using DIR4, DIR5/TIM5, DIR6/TIM6/TIM6e, and GOCO06s is 0.039, 0.038/0.040, 0.036/ 0.038/0.038, and 0.037 m, respectively. The SD of misfits of the Gaussian filtered MDT increases to 0.041 m. The statistics of misfits show that the estimated MDT using GOCO06s and DIR6 perform the best in this study. However, the MDT based on DIR6 only shows 0.3 cm better agreement with the Synthetic_MDT than that based on DIR4. Besides, when comparing with the same release of the GGM, the estimated MDT using direct GGM solution shows better performance than the timewise GGM solution. For example, the SD of misfits of the estimated MDT using TIM6 show~0.2 cm larger than the MDT using DIR6. Estimated errors are deduced based on the error propagation using DIR6, shown in Figure 5, where cs1 is d/o 260. It is seen that most of estimated errors are found near the coast region, especially in the southwestern part of the research area. Statistically, estimated errors (SD) reach 0.1 m near the coast area, while they reduce to less than 0.05 m in most of the open sea.
Further, the effect of different choices of the cut-off frequency for cs1 is investigated, and the SD of misfits between the MDT solution using DIR6 against Synthetic_MDT is shown in Figure 6, where cs1 is selected from d/o 180 to d/o 300 with a step of 10 d/o, corresponding to the expansion of DIR6. The blue dash-line and the orange dash-line in Figure 6 shows the improvement of MDTs using LS-based approach and Gaussian filtering with different choices of cs1. The SD of misfits of estimated MDT solution is always smaller than the Gaussian filtered one, and it decreases gradually with the increase of cs1 and the curve of the SD converges at d/o 260, where the SD of misfits is~0.036 m. Misfits for the Gaussian filtered MDTs rapidly decrease from~0.070 m (d/o 180) to~0.039 m (d/o 300). This is due to the signal contribution from GGM with increasing cut-off frequency. However, the curve of the SD starts to converge after the expansion of d/o 260, and the SD of misfits is improved by less than 1 mm when the expansion of cs1 exceeds d/o 260. The reason for this may be because the SNR of GGM is getting close to 1 with the increase of the expansion, and thus less signal is introduced in the modeling .

Geostrophic Current Validation
Geostrophic current velocities are calculated based on estimated MDT and the Gaussian filtered MDT, where DIR6 in d/o 260 is used in the calculation. The MDT signal may be smoothed out from the coast or current, and into the surrounding ocean due to filtering, causing current signal attenuation. Accordingly, it is necessary to utilize geostrophic current data for validating deduced geostrophic current using estimated MDT and the Gaussian filtered MDT. Geostrophic current velocities can be given as (Lagerloef et al., 1999): where u and v denote the eastward and the northward component of geostrophic current velocity, respectively. f denotes the Coriolis force coefficient, and R, ϕ, λ, and g denote the mean radius of the Earth, the latitude, the longitude, and the acceleration due to gravity, respectively. In the top panel of Figure 7, estimated geostrophic currents deduced from different MDT solutions are illustrated, where First, geostrophic velocity information from Synthetic_Current assimilated a large number of in-situ measurements to enhance the current signal in the model. In addition, the expansion of GGM and the imposed spatial filtering or smoothing condition will make an attenuation of geostrophic current signal. It is notable that the signal of geostrophic current in the open ocean in Figure 7B is relatively larger compared to Figures 7A,C. The main reason is that the MDT estimated using the LS-based approach is still not as smooth as the Gaussian filtered MDT, so the geostrophic current signal is amplified where the gradient of the MDT is overestimated.
It is seen from the bottom panel of Figure 7 that both deduced geostrophic current fields are relatively weaker than the Synthetic_Current. For instance, geostrophic current velocities derived by the Gaussian filtering approach show significant misfits in the coast of the Gulf Stream. However, by comparing Figures 7D,E, it can be seen that the signal attenuation problem in geostrophic current obtained by the LS-based approach is better suppressed, where misfits are smaller. More importantly, the deduced geostrophic current velocities using Gaussian filtered MDT shows stronger noise signal in the coastal area in the southwestern part of the research area, indicating a signal distortion problem. Estimated MDT recovers better geostrophic current velocities in both the coastal and current regions than the Gaussian filtered MDT. Table 6 shows the statistics of misfits between geostrophic current obtained from the LS-based approach and Gaussian filtering, against the Synthetic_Current. It can be seen that the SD of misfits of geostrophic currents derived by the LS-based approach is 0.062 m/s, which shows 1.8 cm/s better agreement than the Gaussian filtered MDT.

CONCLUSION
The traditional geodetic MDT modeling approach requires a filter procedure in the spatial or spectral domain to derive a reasonable solution, which inevitably lead to signal attenuation problems. In this study, we focus on using an LS-based approach for MDT modeling to limit the usage of spatial filtering. This approach constructs observation equations segmentally and carefully models the power matrices of each spectral domain. Further, the unknown parameters of the MDT are estimated in the LS system. Satellite-only GGMs including GRACE-only GGM and seven GOCE-based GGMs with error information are used for MDT modeling, where the effects of the error variancecovariance information from GGM and the contribution of the geoid signal in cs1 are carefully investigated. Finally, deduced geostrophic current velocities based on the MDT solutions are validated and analyzed. Results show that applications of full error variancecovariance matrix or only diagonal error variance of GGMs may have marginal effects on the MDT modeling. The latest GOCE-based GGMs, i.e., DIR6 and GOCO06s, reveal the best geodetic MDT solution using the LS-based approach, where the SD of the MDT misfits is 0.036 and 0.037 m if d/o 260. The estimated error map based on the error propagation illustrates that errors are more concentrated near the coastal region (~0.1 m in Figure 5), while it reduces to~0.05 m ( Figure 5) in the open ocean. The validation using deduced geostrophic velocities shows that estimated MDT recovers better geostrophic current signal in both the coastal and current region than the Gaussian filtered MDT. The SD of geostrophic current misfits by estimated MDT shows an improvement by 1.8 cm/s, comparing with the corresponding ones estimated by the Gaussian filtered MDT.
More future work need to be done to further assess and improve the LS-based MDT modeling approach. First, in order to study the variation of ocean currents due to the global sea level rise over the years, it is necessary to adjust the reference time of the MSS to the last 20 years (e.g., 1997-2018). It is notable that the MDT solutions derived in this study are referenced to the time period 1993-2012. However, with the study of the accumulated altimetry data, it is  found that the sea level is rising gradually due to global warming (Hamlington et al., 2020), which will significantly affect the MDT and deduced geostrophic current. Research reveals that there is more prominent sea state variation near ocean current. For example, Gonçalves Neto et al. (2021) showed that the MSS in the Gulf Stream increased by an average of 4.5 cm from 2009 to 2018 compared to 1993-2007, and there is a significant increase (up to 15 cm) in the northern part of the Gulf Stream and a decrease of up to 15 cm in the southern part of the Gulf Stream. This means that the Gulf Stream is moving to the north. Therefore, the MDT and geostrophic current signal referenced to the last 20 years needs to be further investigated. Second, it is important to use independent in-situ data for robust validation of MDT solutions. For example, the drift data from the Argo program should be reprocessed and used as validation data. And in order to facilitate studies on the vertical datum unification, the performance of the estimated MDT in the coastal region needs to be evaluated using tide gauges. What's more, a comprehensive comparison and analysis of several up-to-date MDT modeling methods, e.g., Andersen et al. (2018a), Karimi et al. (2020), and Wu et al. (2021), together with the LS-based approach need to be performed, and their performance near current and coast region should be further investigated. In this study, the spatial resolution of the MDT is lower than that of DTU15MDT/CNES-CLS18MDT due to the computation limitation of the global SHs. One solution is to replace SHs with a radial basis function (Wu et al., 2018), which could reduce the size of the design matrix and save the computation resources. It is also necessary to evaluate the error information of the altimetric data and introduce altimetric data into the LS system for the MDT modeling.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HS designed this research, carried out experiments and data analysis, visualized the figures and drafted the manuscript. YW, OA, and PK conceived and coordinated this research. XH and YL participated in research coordination. ZZ participated in figure visualization. All the seven authors contributed to discuss results and reviewed the manuscript.