*Article* **Key Factors for Improving the Resolution of Mapped Sea Surface Height from Multi-Satellite Altimeters in the South China Sea**

**Lei Liu 1,2, Xiaoya Zhang 2,\*, Jianfang Fei 2, Zhijin Li 1,3, Wenli Shi 4, Huizan Wang 2, Xingliang Jiang 1, Ze Zhang <sup>2</sup> and Xianyu Lv <sup>2</sup>**


**Abstract:** A satellite altimeter measures sea surface height (SSH) along the nadir track. Multiple satellite altimeters have been in orbit, and the measurements been merged for mapping mesoscale eddies of ~100 km in size in the oceans. The capability of the mapped SSH for resolving mesoscale eddies depends on mapping algorithms. A two-dimensional variational (2DVAR) algorithm was implemented to generate mapped SSH at a grid size of 1/12◦ in the South China Sea. A range of comparisons were performed between the mapped SSH and the commonly used AVISO (Archiving, Validation, and Interpretation of Satellite Oceanographic satellite data) mapped SSH data product at a grid size of 1/8◦ and 1/4◦. The effective resolution, which represents the spatial scale that the data can resolve, was examined. The effective resolution of the mapped SSH using the 2DVAR algorithm is approximately 100 km, while it is 250 km with the 1/8◦ and 1/4◦ AVISO data products. The difference in the effective resolution results from the difference in the background state and thus the background error. The result suggests that the effective resolution of the mapped data could be increased by choosing a background state so that the associated errors could have a smaller decorrelation length scale.

**Keywords:** effective resolution; deviation from background field; merged maps; mesoscale eddies

#### **1. Introduction**

Abundant in situ observations of oceans (e.g., Argo, drifter, and glider) have greatly promoted the progress of oceanography and led to unprecedented advancements in oceanic dynamics [1–5]. However, limited by the high cost of deploying oceanic observing platforms, the density of in situ observations is still insufficient to resolve mesoscale eddies at a scale of ca. 100 km in size over the wide-open ocean [6,7]. Mesoscale eddies have long been recognized as "synoptic systems" in the ocean [8]. With a rapid advancement of remote sensing technology, satellite altimeters, which measure sea surface height (SSH), provide most important measurements for detecting and monitoring mesoscale eddies [9].

Satellite altimeters measure sea surface height (SSH) only at nadir points, producing one-dimensional along-track SSH [10]. To resolve two-dimensional (2D) mesoscale features, it is therefore necessary to merge along-track data from multi-satellite altimeters to produce spatiotemporally continuous maps [11,12]. Currently, the most commonly used gridded altimeter maps are from the Data Unification and Altimeter Combination System (DUACS), formerly known as Archiving, Validation, and Interpretation of Satellite Oceanography (AVISO) data products. It is also called level 4 (L4) product. This L4 product was generated

**Citation:** Liu, L.; Zhang, X.; Fei, J.; Li, Z.; Shi, W.; Wang, H.; Jiang, X.; Zhang, Z.; Lv, X. Key Factors for Improving the Resolution of Mapped Sea Surface Height from Multi-Satellite Altimeters in the South China Sea. *Remote Sens.* **2023**, *15*, 4275. https://doi.org/10.3390/ rs15174275

Academic Editor: Kaoru Ichikawa

Received: 25 April 2023 Revised: 25 August 2023 Accepted: 29 August 2023 Published: 31 August 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

using an optimal interpolation (OI) algorithm. The newly reprocessed delayed-time (DT) DUACS daily global absolute dynamic topography (ADT) maps (DT2018 and DT2021) have a global spatial resolution of 1/4◦ × 1/4◦ (longitude × latitude), providing an effective resolution (ER) of ca. 240 km at middle latitudes [12].

An ER represents the spatial scale that the data can resolve. We emphasize that an ER differs from the grid size. Unlike the grid size, an ER is defined as a space–time scale of the structure that the gridded data can correctly resolve. An OI method always imposes a filtering effect [13,14]. One question is then whether there is room to increase the ER beyond the AVISO data product [1,15].

The SCS is an immense marginal sea in the near-equatorial tropics. Energetic mesoscale eddies occur at a range of sizes. More than 7000 eddies were identified from the satellite altimetry data during the period 1992–2009 [16]. Based on their generation mechanisms, the eddies can be separated into four geographic regions [17]. The Kuroshio invasion is particularly considered in the surrounding area of the Luzon Strait. More than half of the eddies' radii are 100–200 km, <15% have radii exceeding 200 km, and they have a Rayleigh distribution peak of 130 km [18]. Eddies with the largest radii (100–150 km) occur in the central and western SCS. The SCS region is thus selected in this study.

A two-dimensional variational (2DVAR) method is implemented to map ADT data. ADT should be used, rather than sea level anomalies (SLA) for eddy detection [19]. The 2DVAR method was recently introduced for mapping the multi-satellite altimeter data for the East China Sea, the SCS, and the California Current system. The merged data from those areas was shown to have better ER than the AVISO data product [11,12,20].

The 2DVAR and OI methods are based on optimal estimation theory. Both methods use a background state. The associated background error covariance thus plays a key role in the methods. The background error correlation can be characterized by its correlation length scale. The main objective of this study is to illustrate that the background error decorrelation length scale is a factor that dictates the ER of the mapped data. To improve the ER, a proper background state should be carefully chosen. The background error should have a relatively smaller correlation length scale.

The ER in 2DVAR mapped data and AVISO data products are compared. The comparison is made not only with the 1/4◦ global product but also with a customized regional 1/8◦ product in the SCS, particularly provided by DUACS. Another focus is to examine whether the specially customized regional 1/8◦ product has the same performance as other regional AVISO products. Finally, this work demonstrates the robustness of the 2DVAR mapped data, in terms of its higher ER and quality. To comprehensively validate the quality of the merged 2DVAR and 1/4◦ AVISO products, we used two additional mapping datasets (the 1/4◦ AVISO and the HYCOM products), comparing their accuracy and ER. These findings on the HYCOM ER could be meaningful for future data assimilation, as the HYCOM product can distinguish 80 km-scale ocean structures.

#### **2. Data and Methods**

#### *2.1. Datasets*

Major datasets that are used in the analyses are summarized as follows.

First, the 1/4◦ × 1/4◦ AVISO data product. It is also known as the DUACS-DT2021 maps. For convenience, it is called '1/4◦ AVISO'. This is the latest version of 1/4◦ AVISO. In this version, the following improvements were implemented over the previous DUACS-DT2018 version [14]:


Second, The 2021 customized version of the 1/8◦ × 1/8◦ grid-resolution DUACS ADT product ('1/8◦ AVISO').

Table 1 lists the improvements in the 1/8◦ regional European AVISO product compared with 1/4◦ global AVISO product, as well as the regional spatial ER [16,21,22].

**Table 1.** Variance and effective resolution (ER) of regional Europe products. The first column is the variance of the differences between 1/8◦ DT2021 regional Europe product and SARAL-DP/ALtika independent along-track measurements from 2016 to 2019. In parenthesis, variance reduction (in %) compared with the results obtained with the 1/4◦ DT2021 Global product. The second column is the spatial ER of the 1/8◦ DUACS-DT2018 product in the regional European seas. In parenthesis, average resolution over the basin.


The 1/8◦ regional European AVISO products provided substantially better spatiotemporal resolution than the global product. For the current study, 1/8◦ AVISO products were custom-produced for the SCS region for the first time. We hope that they will match the European Sea AVISO products in quality.

Third, the daily averaged reanalysis product of a fully 3D, multivariate, variational ocean data assimilation system (the Navy Coupled Ocean Data Assimilation (NCODA) System (HYCOM-NCODA, hereafter 'HYCOM'; U.S. Naval Oceanographic Office), obtained by averaging the data over eight periods each day. The HYCOM reanalysis product was produced by using the 3D-Var method to assimilate the satellite altimeter SSH, in two alternative ways: (1) assimilating the along-track data of the satellite altimeter directly, and (2) assimilating the merged maps of 2D horizontal analysis of SSH [23].

To evaluate the reconstruction achieved using the different models, we first used independent data, including remote sensing Multi-scale Ultra-high Resolution (MUR) Sea Surface Temperature (SST) data and in situ Global Drifter Program (GDP) drift buoy data [24,25]. In addition, the AVISO along-track L3 data—from the Sentinal-3A (S3A) and Jason-3 (J3) satellites—were used for observation-based ground truthing. Owing to their low observation error and stable operation, S3A and J3 are used primarily to analyze the merged field error and the ER and useful resolution (UR).

The selected period was from 1 April to 1 September 2018, in the SCS, covering 0◦–26◦N and 100◦–130◦E. Five altimeters were in orbit in this region and period, providing sufficiently dense observational coverage and corresponding to the customized 1/8◦ AVISO product provided by DUACS.

#### *2.2. A Two-Dimensional Variational Method*

2DVAR is based on the 2D variational principle. The solution is optimal in the sense of minimum error variance [11,12,20]. The solution is obtained by minimizing the cost function:

$$J(\mathbf{h}) = \frac{1}{2} (\mathbf{h} - \mathbf{h}\_b)^T \mathbf{B}^{-1} (\mathbf{h} - \mathbf{h}\_b) + \frac{1}{2} \sum\_{s=1}^{N} \left( H\_s \mathbf{h} - \mathbf{h}\_s^o \right)^T \mathbf{R}\_s^{-1} (H\_s \mathbf{h} - \mathbf{h}\_s^o) \tag{1}$$

where *h* is the analysis field, *h<sup>b</sup>* and *h<sup>o</sup> <sup>S</sup>* are the background and observation fields of the SSH, respectively; *B* and *Rs* are the error covariance matrices of the background and observation fields, respectively; the subscript '*s*' denotes the time level, and '*N'* is the total number of

time levels of observation. *Hs* is the observation operator that maps the background to the observations. The superscripted '*T*' indicates 'transposing matrix'.

The analytical solution that minimizes Equation (1) can be written as:

$$h^a = h\_b + \left(\mathbf{B}^{-1} + \sum\_{s=1}^{N} \mathbf{H}\_s^T \mathbf{R}\_s^{-1} h\_s^o - H\_s\right)^{-1} \sum\_{\text{Scq}=1}^{N} \mathbf{H}\_s^T \mathbf{R}\_s^{-1} (h\_s^o - H\_s h\_b) \tag{2}$$

The AVISO products are based on the optimal interpolation principle, and the solution can be written as follows [26]:

$$h\_x^a = h\_x^b + \sum\_{r=1}^M \mathbf{C}\_{xr} \left[ \sum\_{q=1}^M A\_{rq}^{-1} \left( h\_q^\sigma - h\_q^b \right) \right] \tag{3}$$

where *h<sup>a</sup> <sup>x</sup>* is the SSH estimate on the regular grid *x*, *<sup>h</sup><sup>b</sup> <sup>x</sup>* is the background field, *h<sup>o</sup> <sup>q</sup>* is the observation field, *Arq* is the covariance matrix between observation points *r* and *q*, and *Cxr* is the covariance matrix between the observation and grid points. The solutions (Equations (2) and (3)) can be shown to be equivalent to each other [15].

One factor that the two methods was adapted differently is the background fields *h<sup>b</sup> x*. The 2DVAR model used the merged map of the previous day, whereas DT2010, DT2014, and DT2018 (AVISO products) used mean dynamic topography (*MDT*), which is the *SSH* over a period *T* [seven years (1993–1999), 20 years (1993–2012), and 25 years (1993–2017)] above the Geoid, as follows [21]:

$$MDT\_T = \langle SSH \rangle\_T - Geoid\tag{4}$$

where *SSH* is the height of the surface above the ellipsoid of the reference and includes the *Geoid*. The ADT field on the initial day of the 2DVAR product was obtained by directly interpolating the observation data along a satellite track without a background field. Oneday merging can eliminate the error differences between different background fields. The first day of the merging algorithm was not used in the evaluation as a product in Section 3. The 2DVAR method used 10 d of along-track data before and after the central day (21 d in total) and considered the temporal evolution error. Therefore, there was temporal correlation in the background field. The different period *T* of the background fields used in the 2DVAR and AVISO models will introduce different background errors:

$$
\varepsilon\_{\mathbf{x}}^{b} = h\_{\mathbf{x}}^{b} - h\_{\mathbf{x}}^{t} \tag{5}
$$

where *h<sup>b</sup> <sup>x</sup>* is the background field and *h<sup>t</sup> <sup>x</sup>* is the unknown actual value at the background state time. During mapping, the covariance (*B* in Equation (2) and *Cxr* in Equation (3)) is directly related to and is constructed using the background error *ε<sup>b</sup> <sup>x</sup>*. Simultaneously, the background error influences the sea level anomaly field. For example, the longer the averaged period, the closer the background field is to the actual value, and the more accurate the inter-annual and climatic-scale signals are in the SLA [26]. However, an increase in the time-averaged range of maps as the background field may not improve the covariance matrix, in turn hindering the improvement in the resolution capability of the merged maps. To demonstrate this and to analyze the influences of the background field, different *B* (*Cxr*) values based on the above two methods were simulated. Although it cannot be obtained directly by observing its error, the estimation of *B* (*Cxr*) can be used to solve the function equation.

In the mapping method, the background error covariance matrix *B* (*Cxr*) and the covariance matrix of observation error (*H<sup>s</sup>* in Equation (2) and *Arq* in Equation (3)) are used as the weight parameters that determine the dissemination of information in the estimation field of the observed data. Assuming that the observation error is spatially uncorrelated, *B* (*Cxr*) determines the propagation weight of the spatial observation data. In data assimilation, the information propagation weight primarily represents the filtering characteristics [18]. *B* (*Cxr*) always can be decomposed as follows:

$$B = \sum \mathbb{C} \sum \tag{6}$$

where *C* is the correlation matrix.

Estimating the matrix *B* (*Cxr*) thus requires the length scale of correlation matrix *C*, which is the core parameter characterizing the covariance matrix *B* (*Cxr*), because the length scale extracts significant information from the matrix [27]. The auto-correlation function was used here to calculate the correlation *C*; the distance at which the correlation dropped to 0.6 was used as the length scale of its correlation.

According to the definition of ε in the forward calculation of SSH, the simulation of *ε* for 2DVAR and AVISO in this study was specified as follows: 2DVAR used the merged maps of the previous day as the background fields, whereas AVISO used the 7-, 20-, and 25-year climate-averaged fields from 1993. The background error was calculated by subtracting the along-track daily data from the background field. For the linear along-track datasets, the covariance of each datasets was calculated using other data within 90 km. The covariance of all the along-track data was then averaged, and the covariance of the Gaussian distribution was obtained via ordinary least squares fitting.

The autocorrelation function of a uniform and isotropic one-dimensional spatial random field is equivalent to the power spectral density (PSD), according to the Wiener–Khinchin theorem [28]. For spatially uniform background errors, the spectral representation of the covariance matrix *B* (*Cxr*) is equivalent to the spectral representation of the correlation matrix *C*. For homogeneous and isotropic background errors, the covariance in the spectral space is the variance of each wave number, namely PSD. The correlation based on the autocorrelation function should have physical (filtering) characteristics similar to the PSD of the background error. The degree of filtering of small-scale signals should have the same variation as the correlation: that is, these values should increase together [29].

The following section presents the relationship between the length scale of the correlation and the energy proportions of different scale errors. The Fourier transform method was used to calculate the energy spectral density which represents the PSD of the background errors. To highlight the differences in the energy spectral density between different models, the spatial distance was divided into five scales, and the portion below the minimum distinguishable scale in all models was omitted when calculating the percentages.

#### **3. Results**

#### *3.1. Signal Proportion of Different Scales in the Background*

The 2DVAR method generated the smaller correlation length scale than the AVISO method (Figure 1). For the AVISO method, the correlation length scales increased with time in years averaged for the background. This indicates that, even when a finer spatial grid is used, selecting a too-long averaged period for the background field prevents the reduction in the correlation length scale.

Around the SCS, slightly more than 25% energies of the identified signal in the background errors of the 2DVAR method were for small mesoscale signals (50 km to 150 km), and 20% energies for signals of 150 to 400 km (Figure 2). In the background errors of the AVISO method, small mesoscale (<150 km) signals had a deficient proportion (2–6%) of energies in the 7-, 20-, and 25-year averaged background. Those signals with a spatial scale >400 km accounted for 75% (7 y) and 88% (20 and 25 y) of the total energy, respectively. The 2DVAR method generated more signals that have smaller scales because it uses a day-to-day estimate, giving less weight to low-frequency longer spatial scales, whereas the AVISO methods exhibit the opposite effect at each length scale.

**Figure 1.** The correlation length scale and its Gaussian-fit distribution of the 2DVAR and AVISO maps with different background time windows.

As the length of the background time window increases from day to year, the proportion of energies with large-scale signals remaining in the background errors increased, and the correlation length scale increases. Therefore, using a long-time-averaged background field in the mapping implementation hinders improvement of the resolution capability of the merged maps.

#### *3.2. Evaluation of Accuracy*

To validate the reconstruction of small mesoscale structure and consistency with the actual signal, different merged maps were evaluated against independent SST and drifter trajectory data and dependent along-track satellite ADT data. Tides significantly influence the coastal area of the SCS. The continental shelf, at an average depth of 200 m, was not included in the analyses in this section.

#### 3.2.1. Remote Sensing Evaluation

Although the spatial patterns of anticyclone eddies are not entirely consistent in terms of SST, or are even orthogonal in some weaker eddies [30], SST of mesoscale eddies is often characterized by a warmer center of the anticyclone, similarly, cyclonic eddies often have colder centers. To a certain extent, SST is responsive to small mesoscale dynamic processes in the upper ocean. In this section, we compared the sea surface geostrophic current, inverted from the ADT data, with the MUR SST anomaly (SSTA) to verify the accuracy of the estimated ADT structure [31].

2DVAR single-day surface flow and temperature corresponded highly in terms of refined structural features (Figure 3a). The other three models show this correspondence only at a larger scale, with poor correspondence in terms of the same refined structure (Figure 3b–d). The geostrophic current of the 1/4◦ AVISO model was weaker than that of the other models. Although the 2DVAR model achieved better reconstruction of mesoscale eddies, it may induce more small-scale noise in the map the same SST.

#### 3.2.2. In Situ Evaluation

The in situ GDP hourly drifting buoy datasets [3] is widely used to study small-scale high-frequency dynamic ocean processes. The GDP measures ocean current velocity at a depth of 15 m using a drogue that pulls a surface float underwater with each passing wave. The drogue creates strain against the lower hull of the surface float, thereby reducing the effect of wind. Here, a drifter path in the middle of the SCS, far from the coast areas, was selected to reduce the influence of tides. A drifter with a curved moving path was compared to evaluate its ability to capture the mesoscale structure in the various merged maps.

On 29 May 2018 (Figure 4), the buoy paths were parallel to the contours of the ADT gradients of the four models. On 30 and 31 May 2018 (Figure 4b,c), the buoy paths were almost parallel to the contours of the ADT gradient of the 2DVAR model and perpendicular to those of the other maps. The result of the other models are mostly inconsistent with the buoy data shown in Figure 4d, for 1 June 2018. The relationship between the buoy path and the height field gradient is similar to that in panels c–e of Figure 4. Comparing the 2DVAR height field with that of the stand-alone buoy reveals that the buoy rotated counterclockwise around a small eddy with a low center at SSH, consistent with the eddy characteristics obtained using the 2DVAR method. The two AVISO models and HYCOM have poor reconstructions of smaller-scale eddies, thus resulting in buoy trajectory crossovers with their height fields. Therefore, 2DVAR exhibits certain advantages over the other models in terms of its ability to present smaller mesoscale eddies in the ocean.

**Figure 3.** *Cont*.

**Figure 3.** MUR sea surface temperature (SST) and the inversion of geostrophic current from merged absolute dynamic topography (ADT): (**a**) 2DVAR, (**b**) HYCOM, (**c**) 1/8◦ AVISO, and (**d**) 1/4◦ AVISO on 7 July 2018.

**Figure 4.** Absolute dynamic topography (ADT) maps of 2DVAR, HYCOM, 1/8◦ AVISO, and 1/4◦ AVISO (sorted by column) with the drifting buoy path. The days are (**a**) 29 May, (**b**) 30 May, (**c**) 31 May, (**d**) 1 June, and (**e**) 2 June 2018 (sorted by row, respectively). Blue line: drift path of the buoy on the focal days; red line, drift path during the two days before and after. The arrow indicates the geostrophic vectors from the middle moment of the focal day.

All of the models exhibited lower geostrophic velocities than the drifting buoys (Figure 5). Although the HYCOM data were the most concentrated in the regression curve, the 2DVAR regression curve illustrates a better fit, with a slope closest to 1. The 2DVAR model also had the lowest root mean square deviation (RMSD) for the entire map.

**Figure 5.** Geostrophic velocity scatter plot and linear regression curve of (**a**) 2DVAR, (**b**) HYCOM, (**c**) 1/8◦ AVISO, and (**d**) 1/4◦ AVISO with the drifting buoy. The slope of the dotted black dot is 1. Equation *y* on the upper left corner represents linear regression curve which is the dotted line colored in red. *R* is the correlation of linear regression coefficients. *C* is the correlation coefficient between the merged map and the drifting buoy data, and the root mean square deviation (RMSD) between the merged map and the drifting buoy data.

#### 3.2.3. Along-Track Satellite Evaluation

To calculate the root mean square error (RMSE), correlations, and linear regressions of the merged maps and thus evaluate their accuracy, the L3 along-track data from the S3A and J3 satellites were used as the hypothetical true values [11]. S3A has smaller orbital spacing and higher data density satellite than J3, allowing it to cover a larger area in the full-field mapping of RMSE. To illustrate the relationship between the 2DVAR maps and those generated using the other models, an index score (*S*) comparing the deviations of each product was calculated [12]:

$$S = 1 - \sigma\_{\text{2DVAR}}^2 / \sigma\_{\text{Others}}^2 \tag{7}$$

where *σ* is the mean square error for each model; 'Others' refers to the models other than 2DVAR; and *S* is the degree of deviation of other models relative to 2DVAR. A positive *S* indicates a higher error relative to that of 2DVAR, while a negative S reflects an error less than that of 2DVAR. Table 2 lists the mean RMSE and *S* for the entire map, reflecting the accuracy of the merged maps.

2DVAR had the lowest RMSE for the center waters of the SCS (at only 1–2 cm), followed by the 1/4◦ AVISO model (at ca. 1 cm higher) (Figure 6a,c,e,g). The 1/8◦ AVISO model had an RMSE of ca. 5–6 cmca. double those obtained in the previous two models (Figure 6d,h). HYCOM exhibited the highest RMSE (>20 cm, relative to the S3A data) in the northern SCS (Figure 6b) and <10 cm relative to the J3 data (Figure 6f). The large RMSE obtained for the Pacific Ocean, Luzon Strait, and the Celebes Sea was due to the tidal influences in coastal areas and high eddy kinetic energy (EKE) transport from the Pacific crossing Luzon Strait to Kuroshio [32]. In summary, the mean RMSE for the entire map (Table 2) was lowest for the 2DVAR model, at ca. 0.01 cm lower than that of the AVISO global 1/4◦ model, and

substantially lower than that of the HYCOM model. *S* varied in the same way as RMSE, being highest for HYCOM and lowest for 1/4◦ AVISO model.

**Table 2.** The accuracy of 4 models. The root mean square error (RMSE, in cm) for entire map compared with along-track S3A and J3 satellite L3 data, and the mean error index score *S* for entire map of each product compared with 2DVAR.


**Figure 6.** Root mean square error (RMSE) on the satellite track of merged absolute dynamic topography (ADT) compared with S3A (**a**–**d**) and J3 (**e**–**h**) satellite along-track data: (**a**,**e**) 2DVAR, (**b**,**f**) HYCOM, (**c**,**g**) 1/8◦ AVISO, and (**d**,**h**) 1/4◦ AVISO. Note that the color scale for each map is differ from each other for the uniform of colors.

Due to the strong influence of the Kuroshio, the zonal flow in the surface (shallower than 500 m) of the Luzon Strait is typically westward [32]. Thus, most of the high-frequency mesoscale eddies generated from the Bashi and Northern Philippines Strait were propagated from the longest channel of eddy propagation in the SCS to the southwest of the northern SCS (from the Luzon Strait to ca. 18◦N, 112◦E), often with a long lifecycle [33]. The high RMSE values southwest of Taiwan island in both of the AVISO models reflect the same distribution characteristic as the longest channel of eddy propagation. This result may occur because the increased presence of eddies in the channel increases the spatiotemporal variability in sea surface conditions, increasing the difficulty in reproducing smaller eddies. For the two AVISO maps, the poor scale-recognition resolution resulted in higher RMSE values in the channel. While 2DVAR model benefits from its good scale-recognition resolution, its RMSE exhibited very low consistency with the propagation channel.

The 2DVAR model exhibited a linear regression correlation coefficient (*R*) > 0.95 and had the highest correlation coefficient (*C*) (ca. 0.98) (Figure 7a,e). The correlation between the HYCOM data and the S3A hypothesized true value was only 0.64 (Figure 7b), almost 0.30 points smaller than that of the other three models. This is highly consistent with the differences in RMSE values.

**Figure 7.** Correlations between the merged maps: (**a**,**e**) 2DVAR, (**b**) HYCOM, (**c**) 1/8◦ AVISO, (**d**) 1/4◦ AVISO, and S3A (**a**–**d**) and J3 (**e**–**h**) satellite along-track data. The ratio of the black dotted line is 1; *C* is the coefficient of correlation; the red dotted line is the linear regression curve; *y* represents its equation; and *R* is the coefficient of correlation for linear regression.

#### *3.3. Evaluation of Effective Resolution*

Figure 8a presents the ADT sequence of the along-track sampling points from the four merged maps and J3 satellite data for a randomly selected date (9 July 2018). The track (Figure 8b) is sorted in the sequence P1–6. Among the four merged maps, 2DVAR exhibited the best fit and best reflected the satellite signal amplitude and small-scale information, especially at P1–2. HYCOM performed the worst, with the largest deviation from satellite data, especially at P5–6. The 1/4◦ AVISO model deviated slightly from the satellite data. The deviation between the 1/8◦ AVISO and 1/4◦ AVISO data waveform was almost regular at every sampling point. An offset occurred because the time scale of correlation and the MDT of the customized 1/8◦ AVISO product had been adjusted.

Each product was a mapping of multiple satellite data, including the J3 data. Part of the merged data approximately 370 and 420 points did not match the ADT variance in J3. The reasons may not only be related to the smaller temporal scales of these processes than that used in 2DVAR or AVISO but also to their smaller spatial scales than the background error correlation coefficient scales.

ER refers to scale-recognition resolution, a parameter that relates the recognition capability of the altimeter product to the dynamic signal of ocean eddies. Scale-recognition resolution is defined as the minimum resolvable spatial scale of the signal in the merged map, and its value indicates the smallest sea-surface eddy that can be distinguished in the signal from the perspective of energy spectral density [13].

Scale-recognition resolution was calculated using the definitions of ER and UR [21] based on the wavenumber PSD of the ADT field, using the observational ADT data as a spatial sequence with distance as the independent variable. ER was obtained using the ratio of the PSD of the noise to the signal, i.e., the noise–signal ratio (NSR). UR was obtained using the ratio of the PSD between the estimated value of the merged maps and the satellite along-track signal, i.e., the signal ratio (SR). At NSR or SR of 0.5, the wavelengths (λ) corresponding to their positions are ER and UR, respectively. By comparing local phase differences, using ER reduces the large systematic error generated by the comparison between the different phases, although it may be more affected by noise. UR reflects a comparison of spatial sequence spectral amplitudes that relatively better represent the

energy magnitude of the signal available within the product data. The average value of each grid point along the track during the study period was recorded to obtain the resolution along the track.

**Figure 8.** The absolute dynamic topography (ADT) sequence (**a**) from merged maps along the alongtrack points of satellite J3 on 9 July. The red and blue dashed lines are the separation position between different tracks. The J3 satellite track in the South China Sea (SCS) for 9 July (**b**); red dotted line). The sequence P1–6 indicates the direction of movement of different parts of the satellite tracks and corresponds to the sequence of the sampling fragments separated by the red and blue dashed lines.

Data for the same day as in Figure 8 (9 July 2018) was selected to display the PSD and scale-recognition resolution results. In the wavelength range of 90–400 km, the PSD value of the mapping field deviation of the 2DVAR product was much lower than that of the 1/8◦ and 1/4◦ AVISO and HYCOM models (Figure 9a). The PSD value approximately 70–200 km wavelength of the mapping field of the 2DVAR product was higher than that of the other models shown in Figure 9b. Excluding the 2DVAR model, the other three models did not differ substantially in PSD. The order of ER and UR of the same product is consistent among all products (Figure 9c,d). Relative to the two AVISO models, the 2DVAR model had ER and UR values ca. 30% lower, whereas HYCOM had slightly higher ER and UR values.

Both ER and UR of the 2DVAR model were mostly between 50 and 150 km, with the higher resolutions near the coast or islands (Figure 10a–d). The HYCOM data ER showed an extremely high value, (>200 km) and occupying a large area, with the UR mostly between 100 and 200 km (Figure 10b,f). For the two AVISO models, ER showed 150–250 km in a large area, with some being ≥250 km. For the two AVISO models, UR was concentrated approximately 200 km in most parts of the study area. For all four models, ER was much higher in the vicinity of the Philippine Islands and their coastal waters than else. Owing to the presence of anomalously fluctuating signals (noise), which could be confused with small-scale signals in the ocean, ER had limited usefulness for resolving eddies.

**Figure 9.** The power spectral density (PSD) that calculated by (**a**) subtraction and (**b**) non-subtraction of the J3 data (except for the gray line in (**a**), it was the origin data of J3 satellite) from the four merged maps (2DVAR, HYCOM, 1/8◦ AVISO, and 1/4◦AVISO). (**c**) Effective resolution (ER) and (**d**) useful resolution (UR) of the four merged maps based on NSR and ER, respectively on 9 July.

By definition, the ER value of an area with a high error increases correspondingly, and the minimum resolved scale of the ocean eddy signal in the data increases. Comparing the error statistics described in Section 3.2.3, ER > UR for most of the regions with large RMSE. For example, for the southwestern waters regions with high EKE, around the Dongsha Islands of Taiwan Island [16], ER < UR for the 1/4◦AVISO model, while for the other three models, ER > UR. Based on these findings, for zones with large RMSE values, UR is more effective than ER for identifying eddies.

Determining the ability to reconstruct eddies involves evaluating the scale of the spatial field and the time scale. Based on a previous work [20], we used the time-frequency domain and spatial wavenumber as reference elements to determine the 2D distribution of the energy spectrum of the height field. The merged maps were decomposed via a 3D Fourier transform. In terms of the period, 2DVAR captured the signals much better than the two AVISO models, within 20 d and 80–200 km (Figure 11a,c,d). The frequency wavelength PSD of HYCOM was large in the wavelength range of 200–1000 km and in the period 0–100 d (Figure 11b). The 2DVAR model exhibited a very low proportion of energy for spatial scale of <80 km (Figure 11a) relative to that in the 80–200 km range, owing to

filtering by its variational merging. Although the AVISO models had an energy proportion below 80 km (Figure 11c,d), the method filtered out more energy at the 80–200 scale.

**Figure 10.** (**a**–**d**) Effective resolution (ER) and (**e**–**h**) useful resolution UR of the four merged models of (**a**,**e**) 2DVAR, (**b**,**f**) HYCOM, (**c**,**g**) 1/8◦ AVISO, and (**d**,**h**) 1/4◦AVISO in the South China Sea (SCS) and nearby waters used S3A satellite data as true values; the darker the color is, the larger the ER or UR value is.

The 2DVAR and AVISO models were used as the hypothetical ground truthing to examine the HYCOM data, and a parameter similar to ER was used to score the frequency wavelength PSD (*FWPSDs*), as follows [20]:

$$FWPSDs = 1 - FWPSD(ADT - ADT\_{true}) / FWPSD(ADT\_{true}) \tag{8}$$

A PSD distribution with a score > 0.5 is considered accurate and reliable. For the 2DVAR models, the HYCOM data were reliable in time–space frequency domain > 80 km, regardless of the time span (Figure 12a). For the two AVISO models, HYCOM was considered reliable only for frequency domains > 150 km, and improved as the time scale increased (Figure 12b,c).

**Figure 12.** The frequency wavelength power spectral density (PSD) values of the three merged maps.

Based on these results, the 2DVAR merged maps have the highest quality and lowest ER among the four models. In an early study [20], 2DVAR merged maps had a higher ER than the 1/4◦ AVISO merged maps for the East China Sea region, especially in the open ocean and to the north of Taiwan Island, where mesoscale eddies are relatively large and long-lasting. Here, for the SCS, 2DVAR merged maps achieved better performance than 1/4◦ AVISO merged maps in reconstructing mesoscale oceanic structure, and their ER was as low as 130 km, or even 110 km around the middle latitudes [12]. For the California Current system, the 2DVAR merged maps could resolve smaller scales than the global 1/4◦ DUACS-DT2018 maps [11]. These findings together demonstrate that the 2DVAR model has advantages in terms of its merged map quality and ER, and without region-dependent.

#### **4. Discussion**

#### *4.1. Signal Composition in Background Field and Associated Error*

For the AVISO models, the background-field MDT was obtained using a multiyearaveraged SSH field. Therefore, the signal contained in the background field is strongly time-smoothed, and the error associated with the background field contained more largescale ocean circulation signals, interannual variation, and seasonal variation. At the same time, the 2DVAR background field did not smooth the signal for any time or region, and induct an evolution error to maintain all signals in the merging processes. Therefore, the background error with the 2DVAR method comprised more small and mesoscale signals than that with the AVISO.

#### *4.2. Filtering Effect of Correlation Coefficient Scale in Variational Method*

The wavenumber energy spectral density of the background error and the characteristics of the background error covariance have similar physical (filtering) characteristics.

The filtering characteristics of the PSD associated with the background error decrease with the spatial scale, owing to the localization and intermittence of small-scale systems, and small-scale errors can only represent a certain proportion of the background error. When the spatial scale of the correlation coefficient is larger, the proportion represented

by small scale error is less. For authentic ocean signals, when the proportion of the signal energy of a certain scale is lower, the error of that scale is more difficult to correct in the merged maps.

Based on these filtering characteristics, and owing to authentic ocean signals of the AVISO merged maps, it is difficult to distinguish small-scale oceanic structure, due to the larger correlation length scale with the AVISO method.

#### *4.3. The Scale of Effective Resolution Compared with Eddy Radius*

The correlation length scale was 10 km for 2DVAR and 100 km for AVISO. The 2DVAR ER was 50–150 km, while that of both AVISO models was 150–250 km.

For the eddies in the SCS, the 2DVAR ER range includes the peak wavelength of the eddy radius Rayleigh distribution, whereas AVISO fails to include the peak wavelength. The 2DVAR merged map can reconstruct mesoscale eddies more accurately than the AVISO maps.

#### *4.4. The Restriction of HYCOM and the Advantages of 2DVAR*

HYCOM reanalysis was produced using data assimilation. Altimetry data is first used to estimate temperature and salinity vertical profiles, employing the ISOP (Improved Synthetic Ocean Profiles) algorithm. The estimated temperature and salinity profiles are then assimilated into HYCOM. This assimilation method may limit the impact of altimetry data in the reanalysis:


It is believed that the long background field time window of AVISO is the key factor causing its resolution to decline. The 2DVAR model has additional technical advantages [23]:


Comparing with HYCOM, 2DVAR does not depend on the above assumptions such as 'noise level', directly merges along-track data, and applies the correlation length scale to supplement noise filtering, thereby retaining small and mesoscale signals. Therefore, the obtained 2DVAR merged map provides higher quality reconstructions than the map obtained using HYCOM reanalysis data.

#### *4.5. Limitations and Future Work*

Although the ER of 2DVAR product has been effectively improved, many small-scale processes still cannot be resolved due to the temporal and spatial scales. To increase the density of observation and acquire more valid information, 2DAVR introduced the evolutionary error in the observation error (*Rs* = *R*m + *R*e, the observation error covariance matrix *R<sup>s</sup>* consists of measurement *R*<sup>m</sup> and evolutionary error covariance matrices *R*e) to address the difference between observation time and mapping time [12]. In addition, the wide-swath Surface Water and Ocean Topography (SWOT) mission was launched on 15 December 2022. As a result, the findings of this study can be extended to resolve small-scale features in maps derived using data from new multi-satellite altimeters, including SWOT data [34]. The multi-scale data merging will be tried to improve further the ER of merged maps in the future [35].
