A Gaussian Function Model of Mesoscale Eddy Temperature Anomalies and Research of Spatial Distribution Characteristics

Abstract

Mesoscale eddies are ubiquitous oceanic phenomena and play an important role in ocean circulation, ocean dynamics, and the transport of material energy. Temperature anomalies are a crucial parameter that reflects the state of mesoscale eddies. This study proposes a Gaussian function model to fit the vertical temperature anomaly (TA) profile to facilitate the analysis of variations, and the principle of the model is based on the fact that each TA profile tends to fluctuate around one or more peaks. The model is extracted and validated using Argo profiles within cyclonic and anticyclonic eddies in the Northwest Pacific Ocean spanning over the period from 2002 to 2021. The validation demonstrates that the model can accurately recover the vertical TA profiles with a limited number of parameters. This makes it suitable for analysing the spatial distribution patterns that require a large sample count. The analysis indicates that eddies with different TA profiles have a spatial aggregation effect in geographic distribution. Eddies with lower extreme temperature anomalies, at depths of 200–300 m, are mainly distributed along two bands on the north side of the Kuroshio Extension (KE) and the North Equatorial Current. Eddies with extreme TAs at the deepest depth (500–600 m) are distributed along the KE.

1. Introduction

Mesoscale eddies are circular water currents that are ubiquitous in the ocean, extending tens to hundreds of kilometres over tens to hundreds of days, and serve as a key bridge in the energy cascade between large-scale and small-scale variability in the ocean [1]. Eddies, which can be classified as cyclonic or anticyclonic depending on their direction of rotation, are important for understanding long-term climate change as they play a major role in the transport of freshwater, heat, dissolved carbon, and other biogeochemical tracers in the ocean [2,3], effectively influencing ocean circulation [4], large-scale water distribution, and biology [5]. Mesoscale eddies play a pivotal role in ocean mixing and drive global water circulation [6] with profound effects on climate and ecosystems [7]. Enhanced mixing induced by mesoscale eddies has been observed in regions of the strong Antarctic Circumpolar Current and the West Boundary Current, such as the Kuroshio Current [6]. The role of eddies in these areas is a long-standing and observationally challenging topic and eddy science has therefore become a central theme of modern oceanography in recent decades. The thermal structure of the eddies is critical for heat transport [2,8] and air–sea interactions [9,10], including changes in the intensity of tropical cyclones [11,12]. In particular, eddy-induced heat transport is critical for maintaining the world’s ocean and climate states [4], and researches have illustrated that eddies can induce heat fluxes by trapping and transporting water with certain temperature in their cores [2,8]. The vertical temperature anomalies are crucial for heat transport estimation and related research analysis.

With the progress of the global Argo project, fast, accurate, and numerous temperature profiles can be accessed worldwide [13]. The temperature anomaly (TA) profiles derived from Argo floats surfaced in eddies have been widely used to study the vertical structure of mesoscale eddies [14,15]. These TA profiles have a complex structure and typically exhibit one or more peaks, with fluctuations in adjacent depth layers. Laxenaire et al. [16] analysed the evolution of a single long-lived Agulhas Ring using an eddy tracking algorithm to uncover its hydrological properties. The track of the single eddy was made possible by the large extent (100 to 400 km diameter) and intense altimetry signal of Agulhas Rings [16]. More studies use the composite method to construct the averaged vertical structures [4,8,15,17,18,19], which are conducted by compositing sufficient Argo hydrographic profiles that are concurrently captured by sea surface eddies in certain regions. As an example, Chaigneau et al. [15] constructed the mean eddy vertical structure of the eastern South Pacific Ocean by compositing the average of 420 (526) profiles acquired by Argo floats that surfaced into cyclonic (anticyclonic) mesoscale eddies. Also, the 3-D structures of composite cyclonic eddies (CEs) and anticyclonic eddies (AEs) were analysed using collocated altimetry sea surface height anomalies (SSHA) and Argo profiles in the Kuroshio Extension (KE) region [19], and analysis showed that at the centre of the composite eddy, a cooling of −2.008 °C was observed at a depth of 360 m for the CE, while for the AE, a warming of 1.788 °C was observed at a depth of 410 m.

Regional ocean climate can be estimated more objectively using individual cluster analysis rather than a fixed region approach. Dong et al. [20] divided the KE region into 9 grids to study the spatial characteristics of the corresponding composite eddies. However, the coexistence of eddies with different vertical structures in the same grid may affect the accuracy of the characterisation of the composite eddies. Maze et al. [21] grouped eight classes of vertically coherent heat patterns in the North Atlantic by using an unsupervised classification Gaussian Mixture Modelling (GMM) method for Argo temperature profiles. Meanwhile, its authors also pointed out that the GMM is computationally inefficient and not satisfactory for large datasets (more than 7000 profiles) [21], which leads to a limitation of the dataset size and the dimensionality of the vertical depth layers. Therefore, efficient methods for the analysis of Argo temperature profiles need to be proposed.

Currently, there is no valid mathematical model for the TA profile itself and no study has focused on automatically extracting multiple fluctuation parameters, including extreme values, depth of extreme values, and steepness in a TA profile. This study proposes a Gaussian function model that can be combined with multiple Gaussians (MG) to fit and extract the fluctuation information. The large dataset of eddy TA profiles can be effectively divided based on the physical characteristics of thermal structures, such as extreme values, depth of extremes, and steepness of the profile, with the characteristics being obtained from model parameters. The MG model presented in this paper is not limited by sample size when dealing with large datasets. It can objectively calculate the characteristics of tens of thousands of Argo profiles and efficiently perform spatial pattern analysis.

Numerous eddies occur in the Northwest Pacific Ocean (NWPO), where a large number of Argo profiles have also been accumulated, and the three-dimensional structures and transports of mesoscale eddies have been comprehensively investigated using a combination of satellite data and Argo profiles [19,20]. It is the integral contribution of eddies over thousands of kilometres in space and years in time that defines their role in the climate system [4]. Therefore, after applying the MG model to the eddy TA profiles in the region, the spatial distribution of mesoscale eddies with different types of TA profiles was analysed.

This paper is organized as follows: Section 2 presents the details of the study area and the data used in this study. Section 3 presents the proposed model to fit the Argo TA eddy profile. The verification results are presented in this section, while the sensitivity analyses and comparative analyses of the parameters associated with the model have also been carried out. The application based on the model parameters, that is spatial distribution variations analysis, is shown in Section 4. Discussions and conclusions are drawn in Section 5 and Section 6.

2. Study Area and Data Processing

2.1. Study Area

For this study, the NWPO was selected as the study area. Specifically, the area between 105–180°E and 5–60°N, as shown in Figure 1. The hydrographic setting in the NWPO is complex. It is one of the most tropical cyclone-prone oceans in the world [11], and is also strongly associated with ocean features such as eddies [22]. There is a mixed water region, where the warm Kuroshio current from the south mixes with the cold Oyashio current from the north, and a recirculation region in the KE [23]. The KE region has a zonal band of high eddy kinetic energy [24], and the short-term formation of the Kuroshio meander appears to be closely related to these eddy activities [25,26].

2.2. Data and Processing

Altimetry allows eddies to be observed by measuring sea surface heights, where currents rotate around local highs and lows due to a geostrophic balance between pressure gradient forces and kurtosis acceleration. In this study, we use the Mesoscale Eddy Trajectory Atlas product (META3.2 DT allsat), which can be accessed from https://www.aviso.altimetry.fr (accessed on 21 September 2023) [27]. More specifically, this refers to the multimission altimetry-derived eddy trajectories in version 3.2 DT (for delayed-time), computed with upstream data from the “all satellites” constellation maps distributed by the Copernicus Marine Environmental Monitoring Service (CMEMS). The product provides anticyclonic and cyclonic eddies detected from the multimission altimetry datasets, with their location, contours, amplitude, radius speed, and associated metadata, which was produced by SSALTO/DUACS and distributed by AVISO+ with support from CNES, in collaboration with IMEDEA.

The atlas was downloaded on 21 September 2023, and covers the period from January 1993 to February 2022. The algorithm used for the product is derived from Mason et al. [28] and further described in Pegliasco et al. [29]. The dataset employed in this study spans a period of 20 years, from 2002 to 2021, which is sufficient to establish and validate the model dedicated to reconstruction of the vertical TA profiles. Similar to the dataset selection in article [30], this paper has chosen mesoscale eddies with tracking trajectories lasting longer than 10 days to ensure more stable data. The dataset has been accordingly divided into the two types, cyclonic eddies (CEs) and anticyclonic eddies (AEs).

Vertical eddy thermal information could be observed with Argo profiles thanks to the Argo project. The Argo project and its data have contributed to the study and development of marine science [31], which provide a large number of vertical CTD (conductivity, temperature, and depth) profiles from the upper 2000 m of the ice-free global ocean and currents from intermediate depths [13]. The delayed mode Argo float profiles were collected and made freely available by the International Argo Program and the national programs that contribute to it (https://www.seanoe.org/data/00311/42182/), accessed on 1 March 2023. The Argo Program is part of the Global Ocean Observing System. Figure 2 shows the Argo array statistics within the study area on a randomly selected day.

All data collected by Argo floats are publicly available in near real-time via the Global Data Assembly Centers (GDAC) in Brest (France) and Monterey (California) after an automated quality control. Additionally, the China Argo Real-Time Data Centre synchronises with the GDAC data server every two days to provide quality-controlled data. In this study, the Argo profile data of the interested area from January 2002 to December 2021 are chosen, and they are downloaded from ftp://data.argo.org.cn/pub/ARGO/, accessed on 1 March 2023. Before using the profile data, it is necessary to control quality to ensure the research results are more reliable. Only profiles with quality markers of ’good’ or ’probably good’ are retained, and additional data screening following Chaigneau et al. [15] is applied to these profiles with rules listed in

Table 1. Additional data screening rules for Argo profiles.

Items	Condition
Minimum observed depth	<10 m
Maximum observed depth	>1000 m
Number of valid points	>40
Depth interval between 0∼100 m	≤ 25 m
Depth interval between 100∼300 m	≤50 m
Depth interval between 300∼500 m	≤75 m
Depth interval between 300∼1000 m	≤100 m

.

Then, based on the altimetry-based eddy identification products, the corresponding Argo profile must be mapped. The decision criterion is whether the position of the profiles lies within the effective boundaries of the simultaneously identified eddies. The process is illustrated in Figure 3. It can be seen that the profiles are extracted in time and space over both the matching Argo data and the mesoscale eddies. The profile data is then vertically interpolated at 5 m intervals using cubic spline function to ensure a smooth and accurate representation of the data.

Furthermore, the temperature anomalies of Argo profiles in eddies are obtained by subtracting the climatological profiles from the Argo temperature measurements. The World Ocean Atlas 2018 (WOA18) [32] is a set of objectively analysed climatological fields of in situ temperature, salinity, dissolved oxygen, and nutrients at standard depths for annual, seasonal, and monthly compositing periods for the World Ocean with 102 depth levels from the surface to 5500 m. The 2018 release updates previous versions of the World Ocean Atlas to include approximately 3 million new oceanographic casts added to the World Ocean Database (WOD) since the previous release, as well as renewed and updated quality control.

For research purposes, the monthly statistical mean temperature on the 1/4° grid for the decadal periods 2005–2017 is used in this study. The corresponding data are available at https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/, accessed on 24 February 2024. Similarly, the monthly climate fields from WOA18 corresponding to the position and month of each Argo profiling float are also interpolated using cubic spline function to have 199 depth levels with 5 m depth steps between 10 and 1000 m. Vertical interpolation was used to simplify the model regression method so that all data had the same vertical axis and to avoid occasional gaps. To make best use of the data, and to avoid interference from surface environmental factors, the shallowest depth level was chosen to be 10 m. Four Argo temperature profiles, falling within the CEs and AEs on a single day, are randomly selected and presented in Figure 4a. These profiles exhibit varying near-sea surface temperatures, with minor temperature fluctuations in the mixed layer and decreasing temperatures following different gradients with increasing depth. Furthermore, the climatic temperature profiles corresponding to the same locations and months of these Argo profiles are shown in Figure 4b. The profiles are smoother than those shown in Figure 4a, and once more, different characteristics can be observed. To demonstrate the temperature anomalies induced by mesoscale eddies, it is necessary to perform a comparison operation of the two.

As a result, the TA profile data are obtained by subtracting the WOA18 climatology data at the same or nearby location and month according to each depth level as (1).

$$T_i^{\prime }=T_i-{TC}_i(i=0,1,2\dots \dots 198)$$

(1)

where T and $TC$ are the interpolated vertical temperature profiles and climatic temperatures, respectively, and i corresponds to the depth index of each layer. There are 199 layers with 5 m as vertical step for range of 10–1000 m.

After processing the above steps, 40,659 cyclonic eddy TA profiles and 44,370 anticyclonic eddy TA profiles were obtained for a period of 20 years from the year 2002 to 2021, and their geographical distributions are as shown in Figure 5, counted in 0.5° × 0.5° grid. It can be seen that the obtained TA profile data cover the study area adequately.

3. Method

3.1. Multiple Gaussian Model

In general, normal AEs induce a positive TA while CEs induce a negative TA. However, there are different types of eddies that cause varying vertical anomalies. For example, cyclonic warm-core eddies (CWEs) and anticyclonic warm-core eddies (AWEs) induce positive and negative temperature anomalies at different depths [26].

Based on these research results, most eddies have warm/cold anomalies in the vertical structure, which consist of more than two peaks/extreme points (in this study, $T_E^{\prime }$ denotes the extreme values of temperature anomaly within the eddy profile). To better fit the TA profiles, we have used the well-established and computationally convenient Gaussian function for fitted models as in (2), which fits the bell-shaped curve very well.

$$f\left(x_n\right)=\sum _{i=1}^M{a}_i{e}^{\frac{-{\left(x_n-u_i\right)}^2}{2{\sigma }_i^2}}(n\mathrm{is}0,1,2\dots 198)$$

(2)

where $u_i$ represents the depth of the extreme value of the ’ith’ Gaussian function. $a_i$ and ${\sigma }_i$ are coefficients representing the amplitude and the standard deviation of the Gaussian curve, respectively. M represents the number of combinations of Gaussian functions, which could be 1, 2, or 3 in this research. Furthermore, for profiles fitted with a single Gaussian function, $a_i$ is equal to ${T_{Ei}}^{\prime }$, which implies the extreme values of the profile within the vertically localised region. The steepness of the Gaussian function is determined by ${\sigma }_i$. This value represents the vertical range of influence of the temperature anomaly. For instance, using the 3$\sigma$ criterion (Lajda criterion), we consider the range of influence of the ’ith’ temperature anomaly fluctuation to be [$u_i-3{\sigma }_i$, $u_i$ + 3${\sigma }_i$]. These parameters are crucial for the thermal structure of mesoscale eddies, especially for the vertical distribution and regulation of heat.

The MG model-fitting procedure developed in this study can automatically extract the best-fit combination solution of Gaussian functions for the TA profiles. The processing steps are summarised in Algorithm 1 and described in detail below.

Algorithm 1 Multiple Gaussian Model Algorithm

Require: Eddy TA profiles data as input $[p_1,p_2,...p_m]$∈$P_m$     $p_i$ = $[t_1^{\prime },t_2^{\prime },\dots t_n^{\prime }]$∈$X_n$, n=199   1: for $p_i$∈$P_m$ do   2: Get $Filtered\_p_i$ via Butterworth low pass filter     with N=3, $W_n$=0.2   3: Get $u_1$,$u_2$,$u_3$     via calculation the relative extrema points of $Filtered\_p_i$   4: Build MG model via $\mathbf{TripleGauss}$     //model function is defined as (2)   5: Get $a_i$∈$(a_{1i},a_{2i},a_{3i})$,${\sigma }_i$∈$({\sigma }_{1i},{\sigma }_{2i},{\sigma }_{3i})$ for $p_i$   6: end for TripleGauss:  $p_i$, $(u_{1i},u_{2i},u_{3i})$   7: Get $pMG\_123i$ (the fitted TA profile via MG model),     $(a_{1i},a_{2i},a_{3i})$,$({\sigma }_{1i},{\sigma }_{2i},{\sigma }_{3i})$   8: Get $R^2\_123i$ //$R^2$ between $pMG\_123i$ and $p_i$   9: Get $(a_{1i},a_{2i})$,$({\sigma }_{1i},{\sigma }_{2i})$,$R^2\_12i$ via DoubleGauss($p_i$,$(u_{1i},u_{2i})$) 10: Get $(a_{1i},a_{3i})$,$({\sigma }_{1i},{\sigma }_{3i})$,$R^2\_13i$ via DoubleGauss($p_i$,$(u_{1i},u_{3i})$) 11: Get $(a_{2i},a_{3i})$,$({\sigma }_{2i},{\sigma }_{3i})$,$R^2\_23i$ via DoubleGauss($p_i$,$(u_{2i},u_{3i})$) 12: Compare $R^2\_123i$,$R^2\_12i$,$R^2\_13i$,$R^2\_23i$,      return with optimum combination solution. DoubleGauss:  $p_i$ and $(u_{1i},u_{2i})$ 13: Get $pMG\_12i$,$(a_{1i},a_{2i})$,$({\sigma }_{1i},{\sigma }_{2i})$ 14: Get $R^2\_12i$ //$R^2$ between $pMG\_12i$ and $p_i$ 15: Get $\left(a_{1i}\right),\left({\sigma }_{1i}\right)$,$R^2\_1i$ via SingleGauss($p_i$,$\left(u_{1i}\right)$) 16: Get $\left(a_{2i}\right),\left({\sigma }_{2i}\right)$,$R^2\_2i$ via SingleGauss($p_i$,$\left(u_{2i}\right)$) 17: Compare $R^2\_12i$,$R^2\_1i$,$R^2\_2i$,      return with optimum combination solution. SingleGaussian:  $p_i$ and $\left(u_{1i}\right)$ 18: Get $pMG\_1i$,$\left(a_{1i}\right),\left({\sigma }_{1i}\right)$ 19: Get $R^2\_1i$ // $R^2$ between $pMG\_1i$ and $p_i$ 20: Return $a_{1i},{\sigma }_{1i}$,${R^2}_{1i}$

(1) Apply a low-pass filter to the TA profiles. A Butterworth low-pass filter with specific parameters (filter order N = 3, cut-off frequency $W_n$ = 0.2) is used to smooth out burrs and jitters in the short distance range while retaining fluctuations in the range of more than 50 m. This step aims to improve the accuracy of automatically extracting extreme points and ensure the correctness.

(2) An algorithm for automatically extracting extreme values of the curve was adopted. This involves extracting both the positive and negative extreme values of the filtered curve, sorting them according to their amplitude, and then selecting the three points with the largest amplitude as expected values of the Gaussian function: $u_1$, $u_2$, and $u_3$. These points represent the depths at which the extreme points are located in the profile.

(3) The Python function ’scipy.optimize.curve_fit’ is used to fit a given TA profile with the model defined in (2) and return the remaining model parameters, including $a_1$, $a_2$, $a_3$, ${\sigma }_1$, ${\sigma }_2$, and ${\sigma }_3$.

(4) In addition, when fitting with triple Gaussian functions, the fit is also compared to the combination of two Gaussian functions formed using any two of the expected depth values, including $(u_1,u_2)$, $(u_1,u_3)$, $(u_2,u_3)$. The combination with the highest coefficient of determination ($R^2$) result is selected. Similarly, when fitting with two Gaussian functions, the algorithm will attempt to fit them separately and return the best solution.

Note that the procedure automatically extracts extreme points during step 1. If fewer than 3 extreme points are obtained, the fitting will switch to a Gaussian function with a lesser number. In this study, the $R^2$ (see (3)) is adopted to measure the degree of fit between the model fitted values and the observed set of values. It normalises the results during processing to compare different models. A value of 1 indicates that there is no difference between the predicted or fitted values and the true values, while a value of 0 indicates that the model is fitted as well as it would be if the fitted value is set directly to the mean of the observed data, and a negative value indicates a less effective model.

$$R^2=1-\frac{{\sum }_{i=1}^n{\left({\widehat{T}}_i^{\prime }-T_i^{\prime }\right)}^2}{{\sum }_{i=1}^n{\left({\overline{T}}_i^{\prime }-T_i^{\prime }\right)}^2}$$

(3)

where ${\widehat{T}}_i^{\prime }$ is the model-fitted value for depth level i, $T_i^{\prime }$ is the observed value for the corresponding level, and ${\overline{T}}_i^{\prime }$ is the mean of the observed values for the profile.

To illustrate the fitting procedure, we take the TA profiles of a CE (Figure 6) and an AE (Figure 7) as examples. In these examples, the MG fitting yields M = 3, which represents a combination of three Gaussian functions.

The TA profile in Figure 6a identifies three extreme points: $P_1$, $P_2$, and $P_3$, all of which have negative values. The corresponding depths are $u_1$ = 30 m, $u_2$ = 155 m, and $u_3$ = 375 m, respectively, with $T_{E1}^{\prime }$ = −4.39 °C, $T_{E2}^{\prime }$ = −4.21 °C, and $T_{E3}^{\prime }$ = −2.95 °C, whereas the three extreme points identified in Figure 7a exhibit different polarities. The $P_1$ at a depth of 25 m has a negative value of −5.25 °C, while $P_2$ and $P_3$ at depths of 555 m and 820 m are both positive, with values of 5.107 °C and 3.13 °C, respectively. It can be seen that the TA profiles of actual mesoscale eddies do not ideally consist of a single peak, but are generally a composite of two or more peaks in a profile curve, where the polarity (positive or negative) of the peaks may be the same (Figure 6a) or opposite (Figure 7a). The three extreme points were found, combined with the low-pass filtered curves firstly, then the parameters of each individual Gaussian function were obtained, and at the end, the parameters were automatically adjusted when fitting via the MG model, which helped to fit the TA profile more accurately.

The extreme points identified in Figure 6 and Figure 7, which represent the centres of fluctuations at different depths in the TA profile, reflect the vertical variations in temperature anomalies caused by the mesoscale eddies. However, it can also be observed that not all fluctuations are significant. For example, the extreme point in Figure 6d is deeper than the other two extreme points. The Gaussian function fitted from this point has a smaller amplitude and the fluctuations are more muted. Consequently, this point does not contribute significantly to the overall profile. A similar situation is observed in Figure 7d, where the identified extreme point does not represent the main components of the TA profile and only contributes to improving the fitting accuracy. Then how to judge the primary fluctuation of eddy temperature anomalies? The fluctuation with the maximum extreme temperature anomaly does not necessarily indicate the primary components of the profiles. As shown in Figure 6, the first Gaussian function with an extreme point of −4.39 °C (maximum temperature anomaly) is located at a depth of 30 m; however, it does not have the largest contribution. Instead, the second Gaussian function with an extreme of −4.21 °C has the largest $R^2$ value with 0.96, indicating the primary fluctuation curve. The same phenomenon occurs in Figure 7. The second Gaussian function with the extreme of 5.07 °C is the primary component instead of the first one with the largest anomaly of −5.25 °C. The maximum extreme values may be caused by drastic changes in the seasonal characteristics of the sea water near the surface. In addition, the burr phenomenon may cause single or multiple anomalies, which can be affected by the external environment during buoy measurements.

Based on the MG model proposed in this study, it is easy to calculate the $R^2$ for each Gaussian function separately, and determine the primary Gaussian bell-shaped curve of the TA profile by searching for the maximum value of $R^2$, that is the single Gaussian curve that can best fit the TA profile curve is considered to be the primary Gaussian function of it. Then, based on this, the corresponding parameters including a (the extreme value), u (the depth at which the extreme value is located), and $\sigma$ (the width of this fluctuation) are obtained, which is the basis for accurate analysis of the physical properties of the vertical structure. The identification of primary fluctuations for further analysis based on the MG model allows for the exclusion of interfering data and a clear analysis of spatial distribution patterns. Based on the MG fitting results, it is shown that approximately 42% of CE TA profiles have a primary function different from the one determined by the maximum extreme value. For AE, this ratio is about 40%. Thus, analysing the features solely based on the TA profile maximum extremes may have implications. This situation will be further explored and compared in Section 4.1.

Table 2. Number results of the profiles fitted by MG model.

TA Profiles	Total Number	Triple Gaussian Fitted	Double Gaussian Fitted	Single Gaussian Fitted
CE	40,659	37,008 (91.0%)	3567 (8.8%)	84 (0.2%)
AE	44,370	39,571 (89.2%)	4677 (10.5%)	122 (0.3%)
Average		90.0%	9.7%	0.2%

shows the results of all eddy TA profiles fitted using the MG model, following the procedure described in previous subsection. The results indicate that the MG model fit better with three Gaussian functions for 90.0% of the eddy TA profiles, compared to the double-Gaussian or single-Gaussian models. The program is designed to compare the effects of model fitting downwards, replacing solutions with a lesser number of Gaussian functions if they fit better. Only a few (about 0.2%) TA profiles were better modelled with a single-Gaussian function, whereas about 9.7% profiles were better modelled with double-Gaussian.

3.2. Model Parameter Setting

How to determine the parameter M? To further assess the possible impact of function model setting, the following sensitivity experiment is carried out. Combinations of two Gaussian functions and four Gaussian functions are applied to the same eddy TA profiles, respectively. That is, a double-Gaussian model and quadruple-Gaussian model are designed for comparison.

The fitting method employed in the procedure is based on non-linear least squares to identify the optimal parameters. However, during the comparison process, it was observed that a significant number of TA profiles were unable to identify the optimal fitting parameters after an upper limit on the number of iterations. This resulted in the fitting failing when it was attempted with quadratic Gaussian functions. This may be attributed to the fact that the quadratic Gaussian function necessitates the estimation of more uncertain parameters, which in turn requires more observed data points. As a result, only 35,849 (88.2%) CE TA profiles were successfully fitted with quadruple-Gaussian functions, while the number for AE was 37,920 (85.5%).

The TA profiles that can be fitted by quadruple-Gaussian functions, are then fitted using both the triple-Gaussian model and the double-Gaussian model. Two measurements of $R^2$ and the root mean square error (RMSE) are calculated for each profile, and

Table 3. Indicator results for different numbers of Gaussian function.

		Double-Gaussian	Triple-Gaussian	Quadruple-Gaussian
Average $R^2$	CE	0.72	0.81	0.86
	AE	−21.82	0.83	0.87
Average $RMSE$ (°C)	CE	0.26	0.21	0.18
	AE	1.11	0.22	0.18

displays the average of the two indicators. $R^2$ has been illustrated in the previous subsection, and RMSE, which indicates the magnitude of the error produced by the model in the fit or prediction, is shown in (4).

$$\begin{array}{cc}\hfill & \hfill RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left({\widehat{T}}_i^{\prime }-T_i^{\prime }\right)}^2}\hfill \end{array}$$

(4)

where ${\widehat{T}}_i^{\prime }$ is the model-fitted temperature anomaly value at the ith depth level of a profile, $T_i^{\prime }$ is the observed value or the interpolated one based on observation for the corresponding depth level, and n is the total number of depth levels.

The results indicate that the double-Gaussian model is not appropriate for fitting all eddy TA profiles with the average $R^2$ less than 0.8 for CEs, and less effective fitting results for AEs. The $R^2$ of most AE TA profiles is less than zero, which may be due to the fact that the automatic curve extreme value extraction procedure fails to capture the primary fluctuation extremes of the curve. Instead, it may have replaced them with other jitter points that had larger anomaly amplitudes but were not the primary fluctuation components.

The triple- and quadruple-Gaussian models had better fits with $R^2$ values greater than 0.80 and RMSE values less than or equal to 0.22 °C, respectively. The quadruple-Gaussian model provides a better fit than the triple model, which is easy to understand, as the combination of more Gaussian functions results in a more accurate fit. It shows an average increase in $R^2$ of 0.05/0.04 (CE/AE) and a corresponding decrease in RMSE of 0.03/0.04 °C. However, this improvement is limited.

Considering that the computational complexity increases with the number of functions, and the previously mentioned problem that four functions would lead to convergence failures for a large number of TA profiles, a maximum of three functions was chosen for this study. The aim of the study is to identify major fluctuations in the TA profiles rather than to capture all anomalies in detail. Therefore, the parameter M = 3 was chosen for the MG model based on the study objectives, the number of model parameters and the computational complexity.

3.3. Model Validation

A series of validations of the MG are carried out. The RMSE and $R^2$, calculated between the observed curve and the model-fitting curve are obtained to quantify the fit and the accuracy of the model. The RMSE, which is the root mean square error between the fitted value and the observed value, is widely used to assess the accuracy of the fitting model.

As shown in Table 3, where the indicators are calculated based on 35,849 (88.2%) TA profiles for CEs and 37,920 (85.5%) profiles for AEs, the average $R^2$ is high enough with 0.81/0.83 respectively. Also, the average RMSE is lower or equal to 0.22 °C, and the automatic fitting model is accurate enough to stand for the TA profiles.

In this study, $R^2$ and RMSE are measured based on the comparison of each vertical profile rather than the lateral comparison of depth layers. According to the calculation based on the observed profiles and the fitted profile results, the distribution is obtained and showed in Figure 8 for CEs and Figure 9 for AEs.

Figure 8 shows the distribution of the two indicators for the fitting results, covering all CE TA profiles from January 2002–December 2021. As presented in figure, the portion with the light green background shows that 90% of the $R^2$ is greater than 0.60, indicating that the vertical profile curves fitted by the MG model are in good agreement with the original observed ones. Furthermore, approximately 66% of the profiles have the metrics of $R^2$ greater than 0.8. The distribution of the RMSE also reflects the good fit of the model, with 57% of the profiles having an RMSE below 0.22 °C and only 10% having an RMSE above 0.342 °C.

Figure 9 illustrates the situation for AEs. About 70% of the profiles exhibit $R^2$ metrics greater than 0.8, with 90% of these having metrics greater than 0.62. Additionally, 52% of the profiles have an RMSE lower than 0.2 °C, while only 10% have an RMSE higher than 0.35 °C.

The aim of the MG model proposed in this paper is to automatically extract extreme values and fit parameter results. However, it is important to note that in some cases, the extracted extreme values may not be the most appropriate results due to factors in the buoy observation process that can result in jittery and anomalous values. The MG model proposed in this study accurately fits around 90% of the profiles and meets the requirements of scientific research. Parameters are automatically extracted, eliminating the need for manual intervention when dealing with large amounts of profile data. This saves time and effort in subsequent analysis. The model employed in this study is significant as it can serve as a methodological reference for downscaling vertical profiles to lower dimensions. Furthermore, it is able to extract parameters from the profiles themselves that correspond to key physical meanings. This enables a more reasonable determination of the primary fluctuation in the vertical profiles, leading to more accurate scientific research and discoveries.

4. Results

4.1. Spatial Distribution of Eddies

The study and analysis of spatial distribution represents a pivotal aspect of mesoscale eddy-related research. Chen and Han [1] analysed the geographic distribution of short-lived and long-lived mesoscale eddies and found that they reside in largely separate geographic zones under different mechanisms. Dong et al. [20] present a detailed spatial distribution of salt and heat transports by eddy movements in the NWPO, which may be useful for the generation mechanism of specific mesoscale eddies. This study presents a novel analysis of the spatial distribution patterns based on the characteristics of TA profiles. A total of 40,659 MG model items for cyclonic eddy TA profiles and 44,370 items for anticyclonic eddy TA profiles from a period of 20 years (January 2002–December 2021) were used to analyse the spatial distribution pattern of eddies. Since the characteristics of eddies with dual/multiple extremes or TA extremes of different polarities are not the focus of this study, parameters in a single-Gaussian function with the largest contribution to the profile curve are chosen as the main physical characteristics, which helps to exclude the interferences and can obviously extract the regular characteristics of the same type of profiles. According to the distribution of depth levels of primary Gaussian fluctuation, we select the relatively densely distributed depth layers for analysis, including the depths of 60–80 m, 100–150 m, 300–400 m, and 500–600 m. And in order to better analyse the spatial evolution characteristics, we also add a depth range of 200–300 m in this study. The distribution of the depth of the extreme point with maximum contribution to the Gaussian fluctuation of TAs for AE/CEs is shown in Figure 10, and areas in light green indicate denser depth level ranges of distribution at local scale. The depth distribution map of temperature anomalies indicates that the mesoscale eddies in the NWPO exhibit a high density of temperature anomalies distributed in the subsurface layer, with a depth of shallower than 200 m.

With sufficient data, the data can be further screened based on the MG model to analyse the geographic distribution pattern more accurately and clearly. First, the depth where the extreme point is located is selected from the Gaussian function model with the largest contribution to the profile curve, and at the same time, the primary Gaussian model is required to have a large enough contribution to the entire vertical TA profile ($R^2>0.5$ is required in this article). Thresholds were also defined for the other 2 parameters: $\left|a_m\right|\ge 0.5°C$, ${\sigma }_m\ge 50$ m, which were used to further filter the profiles with the procedure illustrated in Figure 11.

It is worth noting that, in addition to the widely known CEs with cold centres and AEs with warm centres, ’abnormal’ eddies include CWEs and ACEs have been found in recent studies [26,33,34]. The abnormal eddies are defined according to the surface temperature, including the temperature difference between the eddy interior and the surrounding background, and the temperature anomaly ratio within the eddy. For instance, Sun et al. [26] analysed a case of CWE in Kuroshio–Oyashio Extension (KOE) region, showing positive TA values with a maximum of 1.09 °C at 73.2 m. Meanwhile, ACEs induce a negative temperature anomaly at the ocean’s upper 250 m. In this article, we try to find the ’abnormal’ eddies according to the TA profile, i.e., the profiles with a primary fluctuation curve defined by the abnormal temperature extreme (the positive extreme for CE TA profiles or the negative extreme for AE TA profiles). Stricter screening rules were applied to filter the CWEs and ACEs: $R^2>0.65$ with ${\sigma }_m\ge 50$ m, $a_m\ge$ 1.5 °C for CWEs and $a_m\le -1.5$ °C for ACEs.

Then, the surface location where the Argo profile is located is drawn based, on the region of the depth of the extreme point. Finally, the evolution of the geographical distribution of the different eddy types and the distribution of the ’abnormal’ eddy type are clearly illustrated by plotting the sea surface position of Argo in each 0.5° by 0.5° grid (CEs, CWEs in Figure 12 and AEs, ACEs in Figure 13).

Qiu. et al. [35] have analysed the enhanced mesoscale eddy variability concentrated in two well-defined bands: one along the Kuroshio and its extension paths east and southeast of Japan, along ∼35°N, and the other centred on ∼22°N east of Taiwan, associated with the Subtropical Countercurrent (STCC). As can be seen in the figure, the spatial distribution at different depths of the extreme of the primary Gaussian function eddies also shows a regular variation and pattern.The eddies are relatively sparse, with shallow extreme points in TA profiles (in the depth range of 60–80 m). They have not yet formed concentrated belts, where the northern band is mainly located in the Subarctic Current (SAC) and north of the KE, including the mixed water region, while the southern band is located in the STCC and its southern coast.

It is interesting to note that the distribution of eddies with temperature anomaly extremes in the 100–150 m range shows three bands, with an additional band along the North Equatorial Current (NEC) compared to eddies with depths in the 60–80 m range. With increasing depth, the eddy distribution returns to a two-band structure with increasing concentration, and the southern band moves southwards into the region between the STCC and the NEC. As the depth increases to 300–400 m, the northern distribution shifts to between the area of the Kuroshio Current (KC) part, north of 30°N, and the KE part, south of ∼35°N. It is noteworthy that the eddies with temperature anomaly extremes located between 500–600 m are clearly distributed in the region along the flow axis of the KE.

The distribution of AEs is similar to that of CEs, but when the extreme depth of the maximum contribution function of CEs located at 200–300 m and 300–400 m, their northern zone position is more northward than that of CEs at the same depth. In other words, in range of 200–300 m, their distribution moves northward of the KC and KE, located between 35–40°N, while moving along ∼35°N for CEs. At 300–400 m, their main distribution has reached along the KE, while it is south of the KE for CEs. This is consistent with the pattern of distribution of large amplitude signals detected with altimeters by Itoh and Yasuda [24]. That is, the zones of highest amplitudes are located north and south of the axis of the KE for AEs and CEs, which represent warm-core and cold-core rings, respectively.

ACEs (Figure 13f) are located near the western boundary of the Bering Sea, along the Kamchatka Peninsula and the Japan Trench, which is consistent with the distribution analysis results of Sun et al. [26]. Some are also scattered along the KE and northern part, including mixed waters, and it is noteworthy that a significant number of ACEs are also found along the NEC. The distribution of CWEs (Figure 12f) is similar to that of ACEs, except that it is more concentrated along the KE, and there are also scattered distributions in the South China Sea.

The occurrence and behavior of mesoscale eddies are known to be related to both bottom topography and current variability, including meandering of the KE, southward intrusions and eastward extension of the Oyashio [24]. The interactions and principles between ocean currents and mesoscale eddies have been studied from several perspectives. Both satellite observations and numerical simulation have suggested that short-term Kuroshio meander formation is triggered by anticyclonic eddies that originate in the KE [36]. Qiu and Chen [35] have revealed that baroclinic instability of the vertically sheared STCC−NEC mean flow is the energy source for the regional enhanced eddy signals along the STCC band. The geographical distribution of eddies based on the parameters from the MG model in this study can clearly show the geographical locations and evolution of different types of eddies, which provides supporting material and ideas for the research of the interaction between eddies and ocean currents.

5. Discussion

Using the MG model proposed in this article to fit the TA profiles allows for precise extraction of their main features, such as extreme depth and amplitude. By setting thresholds on the model parameters, the profiles can be classified and filtered for clearer analysis of spatial patterns. Section 3.1 has already concluded that relying solely on the maximum values of the TA profiles is not entirely accurate. Around more than 40% of the profiles exhibit significant fluctuation curves that are not introduced by the maximum extreme point. For comparison, we identified the depths of the extreme points with the highest amplitude and illustrated the distribution in Figure 12 and Figure 13. Both CEs and AEs were analysed for two typical ranges: 60–80 m and 300–400 m. The spatial distribution of the profiles with maximum extreme points at different depth levels is shown in Figure 14.

Compared to the pattern shown in Figure 12a and Figure 13a, Figure 14a,c illustrate that the spatial distribution of the profiles, with the maximum extreme point in the 60–80 m range, is relatively dispersed, with no sign or pattern of banding distribution. Similarly, the distribution of profiles between 200–300 m is not concentrated. In comparison to the results of the MG model analysis in Figure 12d and Figure 13d, which show two clearly separated bands, Figure 14b,d show that there are a considerable number of points distributed between the two bands. In particular, Figure 14b does not exhibit a distinct southern band zone as seen in Figure 12d, and Figure 14d does not demonstrate a clear northern band zone as seen in Figure 13d. Comparing patterns of change at different depths would be challenging.

Ocean dynamics distribute heat in space vertically and horizontally. TA profiles reveal the physical distribution of heat along the vertical axis. Therefore, the analysis of the spatial distribution patterns of the TA profiles would support the scientific analysis of heat transport, the oceanic temperature structure, the dynamics of specific oceanic mesoscale eddies, and climate projections related to vertical heat changes.

6. Conclusions

This study proposes a new fitting model for the vertical temperature anomaly of the ocean eddies using Gaussian functions. The temperature anomalies of mesoscale eddies are extracted by matching satellite remote-sensing altimeter data, Argo profile data, and climatological state temperature profiles. After fitting with the MG model, the vertical characteristics of mesoscale eddies including spatial distribution patterns can be further analysed conveniently and accurately using a large amount of long-term data. Compared with most of the previous studies, which used the features of the sea surface data to analyse the eddies, the internal features of the eddies can be more accurately reflected by using the TA profiles [37].

Based on the MG model and the analysis of the spatial distribution of eddies at different depths of extreme temperature anomalies, it is clear that the eddies are spatially aggregated, geographically separated and related to ocean currents. The eddies at lower depths (60–80 m) of extreme temperature anomalies are mainly distributed along the SAC and the STCC. At depths of 100–150 m, the distribution of eddies increases by a band along the NEC. With increasing depth, the distribution of eddies shows a concentration of two bands, one between the STCC and the NEC, the other along the north/south side of the KE, and at depths of 500–600 m, the eddies are distributed along the KE.

Abnormal eddies, including ACEs and CWEs, can be easily obtained through the MG model parameter screening. They are typically located near the western boundary of the Bering Sea, along the Kamchatka Peninsula and the Japan Trench, with some scattered along the KE and northern region. It is noteworthy that a significant number of ACEs are also found along the NEC, and CWEs are more concentrated in the KE region.

Finally, these results provide some possible applications for the analysis of TA profiles of mesoscale eddies, and the application can be extended to other oceanic regions, including the global ocean beyond the Northwest Pacific.

1. The proposed model fits the Argo TA profile well with fewer parameters, reflecting the physical properties of the profile. The method can be used as a reference for downscaling vertical profiles to lower dimensions, making it suitable for processing large amounts of oceanic data and conducting problem analysis and research.

2. The results may also be useful for classifying or clustering eddies in ocean regions, including the NWPO, according to their vertical temperature structure.

3. Furthermore, the model could be employed to analyse the long-term variations in the eddy TA profiles in a specific area. The presence of various types of profiles within a region may affect the results of the averaging/composite analysis, using the MG model. However, the MG model could be used to select a certain type of profile, which would help to obtain precise temporal patterns.

However, it needs to be noted that the MG model focuses only on the fluctuating features of the TA profile, while a more accurate model may need refinement for the subsurface layer, which has more complex variations, with further reference to sea surface temperature information.