Medium-range forecasting of oceanic eddy trajectory

ABSTRACT

The capability of continuous eddy tracking has been significantly advanced over the past decades, thanks to the availability of sea surface height measurements routinely collected by satellite altimeters. However, due to the complexity of eddy kinematic mechanisms and underlying oceanic processes, a medium-range forecast of individual eddy trajectories is a challenging task. This study develops an algorithm to effectively achieve a precise prediction of eddy propagation with a ten-day lead time by building a physics-embedded temporal evolution model. The theoretical eddy phase speed, current field, and eddy properties are incorporated as physical constraints. Compared to only considering eddy past/climatological velocities, the forecast accuracy of embedded physics constraints is improved by 21%. Forecasting of 2.3 million altimeter-derived eddy tracks globally shows that trajectory predictions are valid for eddies moving towards any direction, and the forecast error achieves 9.9 km on a ten-day average. After systematically evaluating the global forecast error, we conclude that the forecasting accuracy is geographically correlated and affected by external forcings (topography and wind) and eddy properties (geostrophy and life stage). Achieving accurate eddy trajectory prediction is expected to provide critical guidance for in-situ eddy observations, and bring substantial improvements for understanding oceanic substance transportations and biogeochemical processes.

KEYWORDS:

• Mesoscale eddy
• medium-range trajectory forecast
• physics-embedded model
• quasi-geostrophic theory
• altimeter and numerical model products

1. Introduction

The generic term ‘eddy’ refers to a kind of oceanic storm, with rotational and turbulent flows of vortices, swirls, rings, meanders, and some filaments and wakes (Robinson 1983). Oceanic eddies with characteristic spatial scales of tens to hundreds of kilometers, and life spans of weeks to several years while propagating hundreds to thousands of kilometers across the ocean, playing a significant role in transporting climatically important properties such as mass, heat, carbon, and nutrients (e.g. Gaube et al. 2014; McGillicuddy 2016; Zhang, Wang, and Qiu 2014).

Eddy detection based on the sea level anomaly (SLA) became possible in the 1990s with the concurrence of two-satellite altimetry missions (e.g. Le Traon and Ogor 1998). The merged product of the coordinated tandem satellite mission provided high enough data resolution to enable the detection of mesoscale eddies (Chelton and Schlax 2003). The spatial resolution of the merging of two or more satellite products had been interpolated onto a (1/4)° × (1/4)° grid, and the temporal resolution had been increased from a 7-day interval to daily over the global ocean, allowing an oceanic eddy to be individually tracked during its lifetime (Le Traon and Dibarboure 2004; Pascual et al. 2006). Over the past decades, a growing development has been made to automatically detect and track eddies (e.g. the Okubo–Weiss parameter method, Chelton et al. 2007; the winding angle method, Chaigneau, Gizolme, and Grados 2008; the flow direction-based methods, Nencioli et al. 2010; the SLA-based methods, Chelton, Schlax, and Samelson 2011). Compared with other eddy identification methods, Chelton, Schlax, and Samelson (2011) pointed out that the SLA-based method can avoid excess eddy detections and additional noise sources. In addition, some studies have illuminated the potential superiority of manual eddy tracking over automated methodologies. For example, Pilo, Mata, and Azevedo (2015) pointed out that eddies can escape detection in sparsely sampled regions or after the process of merging and splitting. These events can distort the eddy, failing to satisfy the algorithm identification criteria. While manual tracking of eddies has shown promise in capturing eddies for a much longer lifetime. Hughes and Miller (2017) developed a manual eddy tracking method and reported several modon eddy pairs whose propagation speeds are super-fast when geared together. Though the automatic detection of eddies in the SLA field has uncertainties and imperfections due to both the interpolation of the gridded SLA product and the detection algorithm (e.g. Amores et al. 2018; Chen et al. 2021a; Chen, Chen, and Huang 2021b), the SLA-based method is constantly being optimized and has been demonstrated for providing credible enough eddy recognition, and has become the dominant eddy detection and tracking algorithm (e.g. Faghmous et al. 2015; Liu et al. 2016; Mason, Pascual, and McWilliams 2014; Pegliasco et al. 2022).

Given the significant progress of SLA-based eddy tracking, eddy propagation has been individually and statistically analyzed in recent decades and highlighted as having a prominent role in transporting heat and substances across the global ocean (e.g. Chelton, Schlax, and Samelson 2011; Fu 2009; Morrow et al. 2004; Zhang, Wang, and Qiu 2014). Yet the forecast of eddy propagation is still a challenging task. It is of great scientific importance for gaining insight into eddy dynamics and transportation, and of practical significance for designing ocean observing systems, fishing programs, and probing underwater acoustics (Wang et al. 2021). The difficulty in precisely forecasting an eddy trajectory mainly relates to its frequent deviation from the theoretical propagation affected by coupled mechanisms (Chen, Chen, and Cao 2022). In addition, the complex process of splitting and merging that occurs during eddy propagation can also pose a challenge to accurately predict their trajectories. A deeper physical understanding of what controls eddy propagation and integrating it into eddy track prediction will be an important step to solving the problem. As briefly summarized below, several investigations have explored the dominant physical mechanisms that control eddy motion in both zonal and meridional propagations. The theoretical framework proposed by McWilliams and Flierl (1979), and Cushman-Roisin, Chassignet, and Tang (1990) documented that nonlinear mesoscale eddies are expected to propagate westward at approximately the long Rossby wave speed driven by the planetary gradient of the Coriolis effect. The advent of TOPEX/Poseidon (T/P) sea surface height (SSH) data prompted Chelton and Schlax (1996) to employ a Radon transform to characterize phase speeds of SSH and firstly noted that the standard theory for freely propagating linear baroclinic Rossby waves may significantly underestimate the observed eddy phase speeds. In parallel, Killworth, Chelton, and Szoeke (1997) pointed out that this major change would be caused by the presence of baroclinic mean flows. Consequently, endeavors were undertaken to reconcile the theoretical predictions of eddy phase speed with observational estimates by incorporating background flow effects (e.g. de Szoeke and Chelton 1999). Several recent studies based on both observations and theoretical arguments have demonstrated that eddy propagation characteristics can be well-approximated by the long-wave limit of the baroclinic mode relative to the background mean flow (e.g. Klocker and Abernathey 2014; Shi et al. 2018; Tulloch, Marshall, and Smith 2009; Wang, Jansen, and Abernathey 2016). Klocker and Marshall (2014) conducted a comparison between observed zonal eddy speeds and theoretical Rossby wave phase speeds, incorporating the Doppler shift induced by the depth mean flow, revealing a substantial alignment between theoretical predictions and observed velocities. This theoretical prediction of the eddy phase speed serves as a fundamental guide for developing eddy propagation forecasts. In addition to the zonal propagation, it has been recognized that the baroclinic eddy propagation would have a large meridional component due to the rotational advection of the background potential vorticity (PV) (McWilliams and Flierl 1979). In subsequent decades, Vandermeirsch, Morel, and Sutyrin (2001) found that besides this direct effect on eddies, a vertically sheared current associated with a baroclinic PV gradient can induce an additional eddy displacement, which could drastically modify the eddy trajectory. This effect was referred to as the ‘baroclinic β effect’ or the ‘effective β vector’ because of its similarity with the planetary β effect (Samelson 2010; Vandermeirsch, Morel, and Sutyrin 2001). In the case of a zero mean flow, the baroclinic β effect reduces to the β effect, while in the presence of a nonzero mean flow, the PV gradient vector can deviate considerably from the meridional. It is further demonstrated in subsequent studies. For instance, Early, Samelson, and Chelton (2011) pointed out that eddy deflection from the westward could be the effect of vertical shear on the total PV gradient vector from ambient currents. Meanwhile, a recent study found that eddy meridional deflection is more geographically correlated rather than polarity determined, which is actually controlled by the PV gradient (Chen, Chen, and Cao 2022). In addition, other mechanisms exist that can affect the propagation of eddies, such as the bottom topography or wind stress (e.g. Jacob, Chassignet, and Dewar 2002; Wilson 2016). Tulloch et al. (2011) have emphasized that atmospheric wind forcing is an essential and direct driver to generate mesoscale variabilities. Khani, Jansen, and Adcroft (2019) diagnosed the subgrid mesoscale eddy fluxes using a spatial filtering approach and found that there are substantial differences in the characteristics of eddy fluxes between regions of westerly and easterly winds in the presence of topography. These studies have highlighted the important influence of wind fields and topographic effects on both eddy properties and motions. However, these physical constraints are poorly utilized in existing eddy propagating forecast models (e.g. Wang et al. 2021).

Several researchers have employed the predictions of eddy propagation by predicting the SLA field based on numerical/artificial intelligence (AI) ocean models (e.g. Ma et al. 2019; Masina and Pinardi 1994; Robinson et al. 1984; Shriver et al. 2007). For example, Hurlburt et al. (2008) observed daily forecast errors ranging from 22.5 to 37 km for eddies in the Arabian Sea and the Gulf of Oman using the 1/32° NLOM model. Xu, Wei, and Yan (2019) demonstrated that the generation and propagation of two typical anticyclonic eddies in the South China Sea (SCS) can be well predicted based on the HYCOM model when their amplitudes are greater than 8 cm. Ma et al. (2019) developed a convolutional long short-term memory network (Conv-LSTM) for 7-day prediction of the eddy evolution, achieving a mean eddy center error of 11.2 km. Notably, they highlighted that Conv-LSTM network predictions exhibited closer alignment with observations compared to those generated by the HYCOM model. Nonetheless, these eddy forecasts based on projected SLA field represent an indirect approach to eddy trajectory prediction, since the predicted results were daily SLA snapshots and eddy tracking was further needed to obtain the eddy propagating trajectory. Furthermore, the predictability of dynamical models can only be increased by better initialization, sufficient data and resolution, and computing power (Oey et al. 2005). Recently, encouraging progress has been made in eddy trajectory prediction in the regional ocean based on the statistical/deep-learning method. For example, Li et al. (2019) constructed a simple statistical model to conduct a prediction in the northern SCS using eddies’ past positions and surrounding geostrophic current fields. The mean flow advection was taken into account by incorporating the altimeter-derived absolute geostrophic flows for the first time. However, they pointed out that these predictors contributed less to the forecast equations compared to the predictors of the eddy past-day speed and the climatological speed. Based on their prediction model, eddy tracks predicted 1–2 weeks in advance coincide well with the observations, with the overall average forecast error of 27.6 km (1 week) and 42.5 km (2 weeks). Similarly, Wang et al. (2020) forecasted eddy propagation for 1–4 weeks based on the long short-term memory network (LSTM) in the SCS, with zonal and meridional errors of ∼29–47 km and 24–37 km, respectively. Subsequently, Wang et al. (2021) proposed a deep learning framework named MesoGRU for 7-day eddy trajectory prediction in the SCS based on 13 predictors, including the eddy velocity, radius, amplitude, and the SLA imagery information. The smallest prediction errors ranged from 8∼18 km in 1∼7 days.

Despite the certain success of the current eddy trajectory prediction algorithm, the accuracy can be improved. For a 7-day trajectory forecast in the SCS, the present error ranges from ∼18 to ∼28 km. This error represents half of the eddy scale (∼60 km in an ellipse morphology, Han et al. 2021), and more than one pixel in high precision numerical model gridded data (e.g. ∼13.6 km in 1/8° grid, and ∼9.1 km in 1/12° grid), which would introduce uncertainties into subsequent studies. The potential for optimizing eddy trajectory prediction remains substantial, provided that the subsequent limitations and shortcomings can be addressed. First, as previously discussed, the incorporated features in the existing algorithms lack physical mechanisms to support accurate eddy forecasts. It is far from sufficient to characterize eddy motions only by their previous velocities and properties. In-depth physical factors, such as the theoretical eddy zonal phase speed and the baroclinic PV gradient, should be considered to provide essential theoretical guidance for eddy trajectory prediction. Second, an effective prediction mechanism for individual eddy trajectories has not been proposed yet. In general, the trajectory prediction approach is to build a model based on all the historical samples in a specific region, and use this model for the prediction of every trajectory within the area (e.g. Wang et al. 2020; 2021). It is understood that eddy motion is transient, for which models trained on large samples can describe the climatological behavior of an eddy trajectory; however, they inevitably have difficulty depicting the transient state of individual trajectories. More importantly, the westward and eastward propagating eddies are regulated by different principal forcing (e.g. β effect on westward-propagating eddies, and current system on eastward-propagating eddies; Chen, Chen, and Cao 2022; Nof 1981; Peng et al. 2021), and considering all eddy trajectories for the forecast model training will inevitably introduce confusion. An accompanying consequence associated with this is that high accuracy forecasting of a relatively long-time span is challenging to implement. Third, the existing algorithms are limited to specific regional domains (e.g. the SCS, Li et al. 2019), thus lacking comprehensive and systematic analyses of the predictability and influencing factors governing global eddy trajectories.

In this study, we aim to design a new eddy trajectory forecast scheme, named the physics-embedded temporal evolution model (PTEM), to provide a highly precise prediction of eddy trajectories in the global ocean. The mesoscale eddies are often described as the ‘weather’ of the ocean, equivalent to the synoptic variability in the atmosphere (Marie et al. 2017). According to the World Meteorological Organization, periods of the weather forecasting classification can be divided into short-range (0–3 days), medium-range (4–10 days), and long-range (longer than 10 days), respectively. Due to the nonlinearity and instability of the weather time-scale features, the dynamical predictions would be inaccurate in a few days, making reliable predictability limited to less than 10 days even after several decades of progress (e.g. ECMWF 2018; Krishnamurthy 2019). Though studies have been reported to provide forecasts of eddy tracks for the next 30 days (e.g. Li et al. 2019), the results have errors of above 100 km which have exceeded mesoscale scale magnitudes of O(100 km). A qualitative assessment considers that the forecast error for mesoscale eddy tracks should be at least one order of magnitude smaller than its spatial scale, in other words, the forecast accuracy of mesoscale processes should achieve the submesoscale magnitudes (∼10 km). In this study, we focus on the medium-range forecast (i.e. 10 days) eddy trajectory prediction. It is noteworthy that this 10-day trajectory forecast aligns with the minimum repeat cycle of currently operated altimeters (e.g. T/P, Jason-1/2/3), facilitating reliable observation and validation of mesoscale eddies. The medium-range forecast based on the proposed PTEM method is effectively applied to predict over two million altimeter-derived individual trajectories from 1 January 2010 to 31 May 2020. It can avoid the aforementioned limitations as it combines both physical variables related to eddy propagation and individual track characteristics of eddy motions, and allows a comprehensive forecast verification of global eddy trajectories. The rest of the paper is organized as follows. The datasets are described in Section 2. The eddy detection algorithm and the proposed PTEM scheme for eddy trajectory forecast are presented in Section 3. Section 4 implements the algorithm for the medium-range forecast of individual eddy propagation and further validates the influence of the current field on the eddy trajectory forecast. Global verification and precision assessment of successively predicted trajectories are discussed in Section 5. The examination of eddy trajectory forecast ability on the near-real-time SLA product is discussed in Section 6. Finally, conclusions are given in Section 7. Further tests of eddy trajectory forecast derived from the numerical model and near-real-time SLA fields are provided in Appendix A.

2. Data

2.1. Merged altimeter data

The ‘all-satellite’ merged SLA data from multiple satellite altimeters (T/P, ERS-1/2, ENVISAT, JASON-1/2/3, Sentinel-3A, HY-2A, Saral/AltiKa, Cryosat-2, and GFO missions) generated by the Copernicus Marine Environment Monitoring Service (https://resources.marine.copernicus.eu/product/detail/) on a global 1/4°×1/4° grid are used in this study for a 10-yr period from January 2010 to May 2020.

2.2. Current data

Three current datasets ranging from 1 January 2007 to 31 May 2020 are included and checked for forecasting performance: (1) the three-dimensional (0–5500 m) gridded reanalysis current product from the GLORYS12V1 with a spatial resolution of 1/12°×1/12°; (2) the gridded absolute geostrophic current calculated using the geostrophic principle from the absolute dynamic topography (ADT), with a spatial resolution of 1/4°×1/4°, and (3) the subsurface current data at a depth of 15 m observed by drifting buoys. The first two current datasets were generated by the CMEMS (available at https://resources.marine.copernicus.eu), and the drifter-measured current data were released by the Drifter Data Assembly Center at the NOAA (available at ftp://ftp.aoml.noaa.gov/pub/phod/buoydata/). Detailed descriptions can be found in Amanda, Marie, and Matthieu (2019), CMEMS (2021), and Rick and Luca (2019), respectively.

2.3. Wind and topography data

The wind and topography datasets are employed to examine the external influences on trajectory prediction accuracy. Daily gridded wind products were released by Remote Sensing Systems (http://www.remss.com/measurements/wind/) from WindSat (1 January 2010–31 December 2019, version 7.0.1; Meissner and Wentz 2009). The topographical data was the ETOPO1 Arc-Minute Global Relief Model obtained from the NOAA (https://www.ngdc.noaa.gov/mgg/global/), and a detailed description can be found in Amante and Eakins (2009).

2.4. World ocean atlas data

The climatological temperature and salinity gridded dataset derived from the World Ocean Atlas 2018 obtained from the NOAA (available at https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/) is used in this study, which is constructed at 102 standard depth levels upper 5500 m in 1°×1° resolution (Garcia et al. 2019).

3. Methods

3.1. Eddy identification and tracking

Mesoscale eddies are detected on gridded 1/4° delayed-time SLA maps using the algorithm developed by Tian et al. 2020. Individual SLA fields are firstly spatially high-pass filtered to remove coherent signals with wavelength scales larger than 20° of longitude by 10° of latitude using a Gaussian filter (Chelton, Schlax, and Samelson 2011). Note that the SLA fields are segmented into 40 regular blocks with a zonal/meridional spacing of 45°/36° for the parallel eddy detection, which has been proven to significantly reduce the computation complexity compared to previous SLA-based methods at the global scale (Liu et al. 2016). Then, we search in these SLA blocks simultaneously for eddy seeds (i.e. local extreme points in each 5 × 5 grid cell) as candidate eddies. SLA contours are calculated at 0.25 cm intervals for levels of the minimum SLA to the maximum SLA to extract eddy boundaries around the eddy seeds (Liu et al. 2016). The amplitude is defined as the absolute value of the difference between the SLA at the eddy center and the boundary along the eddy boundary. Only eddies with amplitude larger than 0.4 cm, and containing more than five consecutive pixels within the contour are retained (Pegliasco et al. 2022). It has been demonstrated that the 0.4 cm threshold of eddy amplitude is a more reasonable compromise between the time of computation and the size of the detected eddies compared to the noise and the resolution of the altimetry maps. Meanwhile, the effectivity of the 0.4 cm threshold as the amplitude criterion has been proven in the construction of various other eddy datasets (e.g. AVISO 2022).

In the next step, an automated tracking procedure is adopted to determine the trajectory of each eddy from the detected eddy dataset. For each eddy identified on Day 1 (t₁), eddies for the next day (t₂) are searched to find the closest one within a search radius of 0.5°. Note that only eddies whose amplitude (A), area (a), and kinetic energy (E) are between 0.25 and 2.75 times of the t₁ eddy are considered as next-step candidate eddies (Tian et al. 2020). Then, a non-dimensional distance (S) is computed to assess the similarity between the t₁ eddy and t₂ eddy candidates (Penven et al. 2005): (1) $S_{t_1,t_2}=\sqrt{{\left({\displaystyle \frac{\mathrm{\Delta }d}{d_0}}\right)}^2+{\left({\displaystyle \frac{\mathrm{\Delta }a}{a_0}}\right)}^2+{\left({\displaystyle \frac{\mathrm{\Delta }A}{A_0}}\right)}^2+{\left({\displaystyle \frac{\mathrm{\Delta }E}{E_0}}\right)}^2}$(1) where Δd, Δa, ΔA, and ΔE are the variations of the distance, area, amplitude, and eddy kinetic energy (EKE) between eddies at times t₁ and t₂, and d₀= 0.25°, a₀  = π×60² km², A₀ = 2 cm, and E₀= 100 cm²/s² are the characteristic values for the eddy distance, area, amplitude, and EKE, respectively (Mason, Pascual, and McWilliams 2014; Penven et al. 2005). Finally, the eddy pair that minimizes $S_{t_1,t_2}$ is selected to be the same eddy that is tracked from the first step to the second step. Ultimately, t₂ eddies with no t₁ counterparts are considered new ones. The tracking process is executed iteratively until all identified eddies are examined. In addition to the mentioned filters on detected eddies, we only retain eddies with a minimum lifetime of 30 days, and propagating with a speed smaller than ten times of the Rossby phase speed. The tracked eddy trajectories are visualized in Figure 1, and the eddy dataset is publicly available at https://data.casearth.cn/en/sdo/detail/5fa668ad1f4600005e005ba3 .

Figure 1. Altimeter-derived mesoscale eddy trajectories tracked globally from 1993 to 2020. Only eddy trajectories that survived longer than 60 days and propagated westward (eastward) with a minimum of 6° (3°) of longitude are presented, and the westward (eastward) propagating eddies are plotted in red and blue (magenta and light blue). AE and CE are anticyclonic and cyclonic eddies, respectively.

It is worth pointing out that both the SLA gridded products and the eddy detection algorithm have limitations. The spatial (∼80–400 km) and temporal (∼10–35 days) decorrelation scales of the merged SLA product will not allow numerous small eddy signals to be captured, especially in the tropical ocean (Ballarotta et al. 2019; Pujol et al. 2016). It is believed that the accuracy of the eddy trajectory forecast proposed in this study will be further improved when better altimeter measurements become available in the future. Meanwhile, the proposed algorithm has also been validated on eddies tracked based on the model SLA fields (see Appendix A), further demonstrating its stability and robustness, and could serve as a support for eddy studies in ocean models.

3.2. Controlling factor estimates for eddy trajectory forecast

3.2.1. Eddy zonal phase speed predicted by the quasi-geostrophic theory

A physical understanding of what controls eddy propagation is an important step to achieving a credible forecast. Theoretically, eddies will intrinsically drift westward at approximately the long Rossby wave speed (Cushman-Roisin, Chassignet, and Tang 1990; Nof 1981). According to the Introduction to Geophysical Fluid Dynamics, quasi-geostrophic (QG) linear planetary waves in a continuously stratified fluid satisfy a dispersion relation of the form (Cushman-Roisin and Beckers 2011), (2) $\omega =-\beta R_d^2{\displaystyle \frac{k_x}{1+(k_x^2+k_y^2)R_d^2}}$(2) where $\omega$ is frequency, $\beta$ is the latitudinal variation of the Coriolis parameter, $k_x$ and $k_y$ are the zonal and meridional wavenumber components. Considering the eddy scales, the appropriate wavenumber was found to be proportional to the average diameter of tracked nonlinear coherent eddies from the eddy census (Abernathey and Wortham 2015). We define the eddy length scale (L_e) as the average radius of all eddies in each 1° × 1° cell. Note that the eddy radius corresponds to the eddy boundary of the closed SLA contour with the maximum average geostrophic current speed. Assuming the eddy stream function can be described by a Gaussian function with an e-folding scale of $\sqrt{2}$L_e (Chelton, Schlax, and Samelson 2011), the eddy wavenumber is defined as $k_x=k_y={\left(\sqrt{2}{L}_e\right)}^{-1}$. The length scale $R_d$ refers to the first baroclinic Rossby radius of deformation calculated following the method by Chelton, deSzoeke, and Schlax (1998): (3) $R_d={\displaystyle \frac{1}{|f(\theta )|\pi }{\int }_{-H}^{0}N(z)dz(\theta >5^{\circ })}$(3) (4) $R_d={\left({\displaystyle \frac{1}{2\pi \beta (\theta )}{\int }_{-H}^{0}N(z)dz}\right)}^{1/2}(\theta \le 5^{\circ })$(4) where $f$ is the Coriolis parameter, and H is the ocean floor depth. $N(z)=(-(g/{\rho }_0)\times \mathrm{\partial }\rho (z)/\mathrm{\partial }z{)}^{1/2}$ is the buoyancy frequency, $\rho$ is the potential density calculated from the WOA 2018 temperature and salinity data using the Gibbs Seawater Oceanographic toolbox (McDougall et al. 2009), ${\rho }_0$ is a reference density (1025 kg/m³), g is the gravitational acceleration (9.8 m/s²), and z is the vertical coordinate. It should be noted that f gradually decays to almost zero in the equatorial area, where the f-plane is transformed to the β-plane (Lagerloef et al. 1999). Eddy movement speeds at low latitudes are dominant by the β effect, and with increasing latitude, the beta β weakens. Considering the wavenumber of eddies, the zonal phase speed predicted by the QG theory can be expressed as, (5) $C_p=-{\displaystyle \frac{\beta }{L_e^{-2}+R_d^{-2}}}$(5) However, in the presence of the baroclinic flow, the standard theory for freely propagating linear baroclinic waves often underestimates the observed phase speeds (Chelton and Schlax 1996). As de Szoeke and Chelton (1999) further pointed out, the second mode baroclinic structure of the background flow facilitates faster wave speeds than standard theory. They explored this in the three-layer model and found that the planetary-wave phase speed is the sum of the vertically averaged mean flow and the phase speed of Rossby waves. It is also evident by the observed reversed propagation (i.e. eastward) of mesoscale eddies over the Antarctic Circumpolar Current (ACC) and the western boundary current regions where Doppler shifting by the background mean flow is significant (Figure 2(a and b)). The advection by the strong eastward mean flow here exceeds the westward propagation speed of baroclinic Rossby waves. Therefore, the linear relation of Equation (5) is considered insufficient to describe the real zonal motions of mesoscale eddies. Klocker and Marshall (2014) argued that eddy phase speeds over most of the ocean are well approximated by the first baroclinic mode Rossby wave phase speed and Doppler shifted by the depth-averaged mean velocity. Shi et al. (2018) pointed out that eddy phase speed is proportional to the intensity of zonal currents with meridional shear. Also, many works have added this new factor to the equation to reveal the eddy propagation relative to the background mean flow (e.g. Killworth and Blundell 2007; Tulloch, Marshall, and Smith 2009; Wang, Jansen, and Abernathey 2016). Therefore, considering the depth-averaged zonal flow (${\overline{U}}^z$), we further obtain the eddy zonal phase speeds by the inclusion of the Doppler shift of the background mean flow, (6) $C_e={\overline{U}}^z+C_p$(6) where ${\overline{U}}^z$ is computed by vertically integrating the vertical variations of the geostrophic velocity related to horizontal potential density gradients based on thermal wind balance (Karsten and Marshall 2002): (7) ${\overline{U}}^z={\displaystyle \frac{g}{f{\rho }_0}{\int }_{z_1}^{z_2}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }y}}dz}$(7) The global distributions of the altimeter-observed and theoretically-predicted eddy zonal speed ($C_e$) are shown in Figure 2(c and d). It is found that, over most of the ocean, the observed eddy zonal propagation is westward except near the western boundaries and the ACC. The most prominent region for Doppler shifting is found in the core of the ACC at observed eastward speeds of ∼2 cm/s, while the theoretically-predicted speed is slightly larger (∼3 cm/s) given uncertainties in the intrinsic Rossby wave phase speed (Klocker and Marshall 2014). By interpreting eddy phase speed for the dispersive dynamics in the presence of instability, the overall pattern of $C_e$ is quite consistent with observations for both direction and magnitude, which will provide important guidance for forecasting eddy zonal motions.

Figure 2. Global distributions of phase map of eddies (a) and depth-mean current (b); Altimeter-derived (c) and theoretically-predicted (d) zonal speed of eddy propagation; Azimuth angle of the equatorward (positive) and poleward (negative) deflections of the mesoscale eddy (e) and the PV gradient (f).

3.2.2. Eddy meridional deflection driven by the baroclinic instability of mean flow

In addition to the Doppler effect, the vertical baroclinic instability of the mean flow would alter the background PV gradient, changing the propagation of eddies, particularly for their meridional motion (Chen, Chen, and Cao 2022; Samelson 2010; Walker and Pedlosky 2002). The reason is that the background current changes the ambient vorticity gradient so that the eddy obtains an extra self-propagation term for the north–south motion (Early, Samelson, and Chelton 2011; van Leeuwen 2007). Therefore, the baroclinic gradient is considered as a factor to reflect eddy meridional deflection. For linear-waves on the large-scale mean flow, the QG interior PV equation is expressed as (Pedlosky 1987), (8) ${\mathrm{\partial }}_x^2\psi +{\mathrm{\partial }}_y^2\psi +{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }z}\left({\displaystyle \frac{f^2}{N^2}{\displaystyle \frac{\mathrm{\partial }\psi }{\mathrm{\partial }z}}}\right)=0}$(8) where $\psi$ is the stream function for the mean flow. PV and surface buoyancy are both conserved by the advection along the geostrophic flow $U=(u,v)=(-{\psi }_y,{\psi }_x)$ (Lapeyre and Klein 2006). The zonal/meridional components of the mean flow PV gradient ($q_x$, $q_y$), determined from the stratified mean vertical velocity, can be expressed as, (9) $q_x={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }z}({\displaystyle \frac{f^2}{N^2}{\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }z})}}}$(9) (10) $q_y=-{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }z}({\displaystyle \frac{f^2}{N^2}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }z})}}}$(10) The vertical variations of the velocity field are computed based on the thermal wind balance (${\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }z}=-{\displaystyle \frac{g}{f{\rho }_0}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }x},{}_{}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }z}={\displaystyle \frac{g}{f{\rho }_0}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }y}}}}}}}$) based on the grided Argo hydrological dataset. We finally computed the total PV gradient $(Q_x,Q_y)$ by integrating the PV gradient from 2000 m (z₁) to the surface layer (z₂), (11) $Q_x={\int }_{z_1}^{z_2}{q}_x$(11) (12) $Q_y={\int }_{z_1}^{z_2}{q}_y$(12) Then, the two orthogonal components are vectorized to obtain the direction and magnitude of the total PV gradient. To reflect the correlation between the meridional drift of eddy propagation and the PV gradient, the meridional deflection angle of the total PV gradient is computed and shown in Figure 2(f), while Figure 2(e) shows the corresponding eddy meridional deflection angle. A general similarity can be found between the deflection angle of eddies and the PV gradient, implying the primary order of the vertical shear of current fields to eddy meridional deflection. In areas where eddies are deflected (e.g. the North Pacific, the Southeast Pacific, and the North Atlantic, Figure 2(e)), the total PV gradient bears a similarly enhanced deflection behavior with a stronger magnitude. Considering the result of this comparison, the magnitude of the total PV gradient (Q) is incorporated into the prediction of eddy meridional motions. It is worth noting that, besides the PVG of the background current, the eddy-induced PVG calculated based on the Ertel’s potential vorticity also plays an important role in affecting its trajectory since this PV anomaly associated with the advection flow (Wu and Emanuel 1995). The incorporation of eddy PVG into eddy trajectory prediction would improve the forecast accuracy to a certain extent (not shown). However, given that the dominant factor of the theoretical Rossby wave phase speed for eddy propagation is derived from linearizing the quasi-geostrophic potential vorticity equation, the combination of background PVG and eddy theoretical velocities is considered to be promising guidance for eddy trajectory prediction.

3.2.3. Observed eddy velocities, properties, and horizontal current fields

In addition to the specific physics-guided predictors described above, the prediction model considers the specific characteristics of each eddy to provide individual forecasted trajectories. First, eddy propagation speed on the previous day ($[U^p,V^p]$), the averaged velocity during the historical modeled days (see detailed definition of the historical modeled days in Section 3.3, $[\overline{U},\overline{V}]$), and the monthly climatological velocity at the previous location of the eddy core ($[U^c,V^c]$) in the zonal and meridional components are included. The speed of zonal and meridional eddy propagation is defined as the ratio of the distance traveled by the eddy over the time elapsed of 7 days, so as to reduce the sampling and accidental uncertainties (e.g. errors associated with eddy centroid identification or eddy-eddy interaction). The climatological velocity is obtained by a vector composition of eddy velocity within each 1°×1° cell, which is expected to minimize the effect of high frequency perturbations induced by instantaneous shift of individual eddies. Second, individual characteristics of eddies including amplitude (A) and eddy geostrophic ratio (G) are included. Eddy amplitude has been found to significantly affect eddy propagation for both zonal speeds and meridional deflection based on simulations and observations. Eddies with larger amplitudes tend to have faster phase speed (e.g. Chelton, Schlax, and Samelson 2011; Early, Samelson, and Chelton 2011; Shi et al. 2018), and show an increasing poleward deflection for both polarities (Chen, Chen, and Cao 2022). Early, Samelson, and Chelton (2011) speculate that the ‘wave drag’ caused by the Rossby waves may be a reason that slows the propagating speed of eddy with smaller amplitude. The eddy geostrophic ratio is defined as the relative length scale of the ratio of the eddy radius and the Rossby radius of deformation, which reflects the deviation magnitude of the eddy scale relative to the Rossby scale. It is an important quantity in determining horizontal scales which characterize eddy sizes and trapping scales, indicating the scale of eddy instability (Stammer 1997). When the eddy geostrophic ratio is less than 1, eddies are more ageostrophic and easily undergo submesoscale instabilities (Chen and Chen 2023). It was found that eddies with a larger length scale propagated slightly faster with a more consistent and stable speed (Chelton, Schlax, and Samelson 2011; Early, Samelson, and Chelton 2011). Our recent work found that the meridional deflection of eddies has a certain correlation with their geostrophic ratio, as eddies with weak geostrophy tend to move poleward irrespective of their polarity. In addition, the spatiotemporal-matched subsurface current field ($[U_F,V_F]$) derived from the GLORYS12V1 reanalysis dataset is introduced. It should be noted that instead of using the current field at a fixed horizontal layer, the depth of the vertical eddy core (i.e. the depth corresponding to the most significant vertical potential density anomalies induced by individual eddies, see Chen, Chen, and Huang [2021b] for detailed procedures of eddy core depth determination) is considered as the horizontal layer with the most significant impact on eddy dynamics.

3.3. Physics-embedded temporal evolution model for eddy trajectory forecast

In this section, we introduce the proposed eddy trajectory prediction algorithm for individual eddy tracks based on the predictors described in Section 3.2, which consists of the following three steps.

3.3.1. Optimizing window size for individual eddy tracks

To ensure better performance of the algorithm, a flexible moving window strategy is utilized to extract the most appropriate individual historical trajectories to be used for fitting the trajectory prediction regression model. To be specific, unlike the existing studies that arbitrarily use a fixed modeling period to predict the entire eddy trajectory sets, we find the appropriate window for each eddy track to minimize the predicted trajectory deviation. As an example in Figure 3(A), excessively long modeling periods tend to incorporate an overload of historical trajectories, confusing the accurate prediction of both the direction and the velocity (Figures 3(A) [a–b]). In contrast, fewer modeling periods may introduce errors due to insufficient samples (Figure 3(A) [d]). To ensure an appropriate number of samples and avoid introducing excessively long historical tracks, the flexible modeling window size is set to a range of 5–15 days. Assuming that the prediction starts on day t + 1, the modeling period ranges from t-n to t-2 (6 ≤ n ≤ 16, orange points), and days from t-1 to t (pink points) are used as the test stages. Let $x_i=[U^p,\overline{U},U^c,C_e,U_F,A,G{]}_i$ and $y_i=[V^p,\overline{V},V^c,Q,V_F,A,G{]}_i$ be the zonal (x) and meridional (y) predictors on day i, the predicted zonal and meridional location, $p^x$ and $p^y$, in the test stages can be expressed as, (13) $[p_{t-1}^x,p_t^y]=f^x(x_{t-n},\cdots x_{t-3},x_{t-2})$(13) (14) $[p_{t-1}^y,p_t^y]=f^y(y_{t-n},\cdots ,y_{t-3},y_{t-2})$(14) where $f^x$ and $f^y$ represent the zonal and meridional prediction regression function described in Equation (20), t-n is the starting point of the modeling window (i.e. n is ranged from 6 to 16 and is iterated in a one-day step, see Figure 3(A)). Suppose $T^x$ and $T^y$ is respectively the altimeter-observed eddy zonal and meridional location in the test stages, the testing error corresponding to each window size is estimated by the great circle distance (GCD), expressed as (Baskar 2021), (15) $GCD=R_e\cdot \arccos [\sin T^y\sin p^y+\cos T^y\cos p^y\cos (T^x-p^x)]$(15) where $R_e$ is the radius of the Earth. When the GCD reaches its minimum, the window size (m) is taken as the fitted window to be the typical trajectory prediction component (Figure 3(A) [c]).

Figure 3. Schematic illustration of the proposed eddy trajectory forecast algorithm. (A) The process of determining an optimizing window size for individual eddy track forecasts. The grey and orange dotted lines represent the historical and modeled eddy trajectories, respectively. The blue and magenta dotted lines represent the two-day forecasted and true eddy trajectories using the corresponding window size. (B) Evolving forecast for medium-range eddy trajectory with ten-day lead time from day t + 1 (e) to day t + 10 (h). The grey and orange dotted lines are the same in (A), the purple lines are the true eddy trajectories and the yellow lines are the forecasted eddy trajectories used as modeled tracks to forecast the next-day eddy location. The green dotted lines are the forecasted eddy position of the corresponding days. (C) The result of the forecasted eddy trajectory.

3.3.2. Evolving forecast for medium-range eddy trajectory

After the modeling window size is determined, the future eddy trajectory is predicted based on the quadratic polynomial regression. Instead of employing the same functional relationship for predicting the future multi-day eddy paths, the function built on a specific historical trajectory only predicts the next eddy stage, and then the predicted eddy stage will be involved in the next model training, so as to develop a temporal evolutionary forecast (Figure 3(B)). Taking advantage of the optimized modeling window size (m), an eddy trajectory can be temporally evolving forecasted for a medium range with ten days lead time from t + 1 to t + d (d = 10). For clarity, we take $C_z=\{(x_z,y_z)|x_z\in \{x_{t-m},\dots ,x_t,x_{t+1},\dots ,x_{t+(d-1)}\},y_z\in \{y_{t-m},\dots ,y_t,y_{t+1},\dots ,y_{t+(d-1)}\}\}$ to represent the feature pair of the zonal and meridional components on day z. In practice, when it is required to introduce modeling features in the future moments (i.e.$t+1\le z\le t+9$), the $C_z$ is obtained based on the previous predicted eddy position and the monthly means of climatological fields. For the sake of distinction, we note this part of features as ${\overline{C}}_z(t+1\le z\le t+9)$, and the predicted eddy location from day t + 1 to day t + 10 can be obtained by, (16) $\mathrm{Day}t+1:P_{t+1}=F_{t+1}(C_{t-m},\cdots ,C_{t-1},C_t)$(16) (17) $\mathrm{Day}t+2:P_{t+2}=F_{t+2}(C_{t-m-1},\cdots ,C_{t-1},C_t,{\overline{C}}_{t+1})$(17)            …  …  (18) $\mathrm{Day}t+9:P_{t+9}=F_{t+9}(\cdots ,C_{t-1},C_t,{\overline{C}}_{t+1},\cdots ,{\overline{C}}_{t+8})$(18) (19) $\mathrm{Day}t+10:P_{t+10}=F_{t+10}(\cdots ,C_{t-1},C_t,{\overline{C}}_{t+1},\cdots {\overline{C}}_{t+8},{\overline{C}}_{t+9})$(19) where $P_{t+j}=[p_{{}_{t+j}}^x,p_{{}_{t+j}}^y]$ is the forecasted eddy zonal and meridional location with j days lead time. $F_{t+j}=[f_{t+j}^x,f_{t+j}^y]$is the zonal and meridional prediction model with j days lead time obtained based on the multiple quadratic regression defined by summing the multiplication between the variables to highlight their characteristics and dynamic interactions. The regression in zonal and meridional are modeled separately for each day. Taking the zonal prediction model on day t + 1 ($f_{t+1}^x$) as an example, the modeling process can be described as follows. Let $\overline{X}=[x_{t-m},\cdots ,x_{t-1},x_t]$ be the zonal predictors for m day on day t + 1, $f_{t+1}^x$ can be summarized by a mathematical formula, (20) $f_{t+1}^x=XW+b=w^1{x}_1+w^2{x}_2+\cdots +w^8{x}_1^2+w^9{x}_1{x}_2+\cdots +w^q{x}_{n-1}{x}_n+w^{q(q+1)}{x}_n^2+b$(20) where $W=[w^1,w^2,\cdots ,w^{q(q+1)}]$ is the weight coefficients, q = 7 is the number of the predictors, $\mathcal{X}$ is a matrix with size m × 56, which contains the primary, quadratic and cross-multiplicative terms of each predictor in $\overline{X}$, and $b$ is the bias coefficient. It's important to highlight that an investigation into the efficacy of employing both primary and cubic functions has been conducted. The results reveal a notable reduction of 16% in the prediction error subsequent to the incorporation of the secondary nonlinear term, indicating the effectiveness of nonlinear quadratic functions in forecasting nonlinear eddy motions. Conversely, the introduction of the more intricate triple term yields negligible alterations in prediction accuracy. Then, the least-squares optimization is used to evaluate the $W$ and b. The main idea is to minimize the squared errors between the true and predicted positions of the tracks within the modeling window. After the optimal coefficients have been determined, the eddy zonal location ($p_{t+1}^x$) on day t + 1 can be obtained by substituting the features of day i + 1 ($x_{t+1}$) into the model $f_{t+1}^x$: (21) $p_{t+1}^x=[w^1,w^2,\cdots ,w^7]{\left[\begin{array}{c}{U}^p\\ \overline{U}\\ ⋮\\ G\end{array}\right]}_{t+1}+[w^8,w^9,\cdots ,w^{49}]{\left[\begin{array}{c}{U}^p\overline{U}\\ \overline{U}{U}^C\\ ⋮\\ G{U}^p\end{array}\right]}_{t+1}+[w^{50},w^{51},\cdots ,w^{56}]{\left[\begin{array}{c}{(U^p)}^2\\ {(\overline{U})}^2\\ ⋮\\ {(G)}^2\end{array}\right]}_{t+1}+b$(21) Likewise, the eddy meridional position on day t + 1 can be obtained from the model $f_{t+1}^y$ based on the meridional predictor y_i. Iterating the above process, the zonal and meridional position of the eddy during time t + 2 to t + 10 can be determined (Figure 3(C)). Figure 4 presents two examples of the forecast of eddy trajectories alongside their corresponding PTEM functions. It can be clearly revealed that for each day of an individual eddy trajectory, the prediction is made by fitting the specific optimal parameters for a function from its previous adjacent historical trajectories. Consequently, the derived predictive values for the given day are seamlessly integrated into the preexisting dataset, thus engendering the formulation of the ensuing day’s predictive function. It is noteworthy to mention that the predictive trajectories in the zonal and meridional dimensions are discretely segregated, employing distinct functional representations. Admittedly, the evolving model process will increase the cost of computing. However, we consider it necessary and effective to improve the trajectory prediction accuracy, especially when eddies deflect in direction or accelerate/decelerate during the predicted path, which is a frequent occurrence in eddy motion under complex oceanic processes.

Figure 4. Example of the medium-range eddy trajectory forecast for an westward propagating AE (a) and eastward propagating CE (b). The bottom panels are the corresponding PTEM model. The $P_{t+j}=[p_{{}_{t+j}}^x,p_{{}_{t+j}}^y]$ is the forecasted eddy zonal (x) and meridional (y) location with j days lead time based on the function $f_{t+j}^x$ and $f_{t+j}^y$. The coefficients in the prediction function are obtained by least squares fitting using the corresponding historical trajectories. The definitions of each variable in the function are detailed in Section 3.2.

3.3.3. Model evaluation

The performance of the proposed trajectory prediction algorithm is evaluated based on eddy tracks in 2018 from two perspectives: the integration of diverse physics-guided variables and a comparative analysis against two recently employed deep learning models, namely, Deep Neural Network (DNN) and LSTM (LeCun, Bengio, and Hinton 2015), as previously applied to eddy trajectories forecast (Wang et al. 2020). Note that the LSTM is built in the classic way of time-series prediction in deep learning (Zhang, Zheng, and Liu 2021), while the construction of the DNN model takes inspiration from the temporal evolution prediction scheme proposed in this study. Three measures are computed to provide a broader view of the forecast accuracy: the adjusted determination coefficient ($R_{Adj}^2$), the averaged great circle distance ($\overline{GCD}$) as defined in Equation (15), and the Pearson correlation coefficient between the true and predicted eddy trajectory. Considering that it is to demonstrate the effectiveness of incorporating physical constraints on track prediction, the $R_{Adj}^2$ is employed instead of the commonly used determination coefficient (R²), given that $R_{Adj}^2$ can effectively avoid obtaining a better model result when one includes more predictors even if the extra predictors are random noise by taking into account the sample and independent feature number ($R_{Adj}^2=1-[(1-R^2)(n_0-1)]/(n_0-p_0-1)$, where n₀ and p₀ are the sample number and feature number, respectively). While the $\overline{GCD}$ and the correlation coefficient are used to assess the prediction performance in terms of the absolute distance error and the morphological similarity between the actual and predicted trajectories, respectively. As summarized in Table 1, the inclusion of physical variables progressively enhances prediction performance. The most significant improvement occurs with the incorporation of physical factors (the $R_{Adj}^2$ increase from 0.76 to 0.88), followed by the inclusion of the current fields (the $R_{Adj}^2$ increase from 0.76 to 0.83). A comparison between Case 0 and Case 6 suggests an improvement of approximately 7 km in prediction accuracy (the $\overline{GCD}$ is reduced from 17.00 km to 10.03 km). It is also evident that the PTEM provides better performance when compared with Case 7 and Case 8 predicted using DNN and LSTM based on the same predictors. Moreover, the PTEM-predicted trajectories exhibit the highest correlation with the actual paths (0.83), illustrating the high consistency between predicted and true trajectory patterns.

Table 1. Performance of the medium-range eddy trajectory forecasting under different conditions.

Case	Forecasting method	$R_{Adj}^2$	Error (km)	Correlation coefficient
Case 0	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$])	0.77	17.00	0.66
Case 1	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [A,A])	0.81	15.72	0.71
Case 2	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [G,G])	0.80	16.02	0.71
Case 3	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$U_F$,$V_F$])	0.83	13.01	0.75
Case 4	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$C_e$,$Q$])	0.85	11.67	0.81
Case 5	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$C_e$,$Q$], [$U_F$,$V_F$])	0.91	11.48	0.83
Case 6	PTEM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$C_e$,$Q$], [$U_F$,$V_F$], [A,A], [G,G])	0.93	10.03	0.83
Case 7	DNN ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$C_e$,$Q$], [$U_F$,$V_F$], [A,A], [G,G])	0.78	15.02	0.74
Case 8	LSTM ([$U^p$,$V^p$], [$\overline{U}$,$\overline{V}$], [$U^c$,$V^c$], [$C_e$,$Q$], [$U_F$,$V_F$], [A,A], [G,G])	0.91	22.11	0.65

Then, histograms of the $\overline{GCD}$ from day 1 to day 10 in different cases are given in Figure 5(a). As expected, the PTEM can reach the smallest prediction error ranging from 3.25 to 19.93 km in 1–10 days (grey lines), while the largest error occurs when using the LSTM method, ranging from 6.67 to 37.38 km (purple line in Figure 5(a)), which further demonstrates the performance enhancement of the proposed algorithm. In addition to the absolute distance error, the comparison of statistically-averaged trajectories is an important aspect to assess the model performance. As shown in Figure 5(b and c), the true and predicted average trajectories of anticyclonic eddies (AEs) in the Northeast Pacific (NEP, Figure 5(b)) and cyclonic eddies (CEs) in the ACC (Figure 5(c)) are calculated. These two regions are chosen to reflect both the westward and eastward propagating eddies with equatorward and poleward deflections, respectively. It can be clearly found that the average predicted trajectories in case 6 almost overlap with the true trajectories for both equatorially-deflected AE and polar-deflected CE. In contrast, the average predicted trajectories significantly deviate from the true trajectory when no physical parameters are taken into account. The necessity and validity of the inclusion of the physical predictors are further demonstrated.

Figure 5. (a) Comparisons of the daily-mean distance errors in different cases; (b) Average of the normalized eddy trajectories relative to the initial location of each eddy for AEs in the Northeast Pacific for the true trajectories, forecasted trajectories in Case 0, and forecasted trajectories in Case 6. (c) as in (b), but for CEs in the ACC. Equator. = Equatorward, Pole. = Poleward.

4. Medium-range forecast of individual eddy propagation

To give an intuitive impression of the eddy propagation forecast result, the medium-range predicted tracks in the Northwest Pacific Ocean (NWP) and the Southeast Indian Ocean (SEI) with predominantly westward and eastward eddy motions are shown in Figure 6(a and b). It is visible that the forecasted eddy propagation is consistent with the altimeter-observed eddy paths for both westward and eastward propagating AEs and CEs. Quantitatively, the ten days averaged GCD between the forecasted and true eddy location in the NWP and the SEI are 5.3 and 6.6 km, respectively, exhibiting slightly lower forecast accuracy in the SEI dominated by eastward propagating eddies. However, the prediction accuracy is not appreciably different for AEs and CEs (5.5 and 6.9 km for AEs, and 5.3 and 6.2 km for CEs in the NWP and the SEI respectively), indicating that the proposed algorithm has no obvious polarity preference. In addition, the scatters between the actual and predicted eddy locations of the longitudes and latitudes appear to fit well (Figure 6(c–f)), suggesting that the proposed PTEM serves as an effective approach in eddy propagation forecasting.

Figure 6. (a–b) Prediction examples of the medium-range forecasts for eddy tracks in (a) the Northwest Pacific and (b) the Southern Indian Ocean. (c–f) Scatter plots of the true and predicted longitudes (yellow scatters) and latitudes (pink scatters) in the Northwest Pacific (c–d) and the Southern Indian Ocean (e–f).

It should be pointed out that the current field, which modulates eddy motions both through the vertical shear and horizontal advection, is the dominant external forcing in regulating eddy movements (Chen, Chen, and Cao 2022), therefore is a critical factor for eddy trajectory forecasting. To provide an in-depth insight into how the background current included in the prediction model affects the eddy track forecast, four typical cases are given in Figure 7. Figure 7(a–d) show the predicted trajectories before and after considering the background current for a westward AE and an eastward CE in the NEP and the ACC, respectively. When the factors of background current are absent, the predicted eddy tracks inevitably tend to follow their historical trajectory for both westward and eastward propagating eddies (Figure 7(a and c)). The predicted trajectories are significantly corrected when the background PV field is incorporated, in which the predicted trajectory almost overlaps with the corresponding true trajectory and moves roughly along the PV gradient contours (Figure 7(b and d)). Meanwhile, quantified error evaluations show that the mean-flow correction reduces the predicted distance error by more than double (23.42 km to 11.52 km for the AE, and 13.10 km to 5.61 km for the CE). Furthermore, the simultaneous horizontal current advection plays an important role in eddy motions. Though the surface geostrophic current derived from the altimeter has been incorporated into eddy trajectory predictions (Li et al. 2019), it is insufficient to characterize the eddy-current field features since eddy motion is accompanied by notable submesoscale, ageostrophic energies (Cheng et al. 2020; Zhang et al. 2019). We approach the problem by introducing the high resolution, subsurface current field in eddy trajectory forecasting: on the one hand, instead of using the geostrophic current with a spatial resolution of 1/4°, a high precision reanalysis current field with a spatial resolution of 1/12° derived from GLORYS12V1 is employed to provide a better approximation to the real flow field. On the other hand, instead of using the surface current field, the subsurface current at the eddy core depth is used considering the current mainly influences the eddy motions at the depth of their vertical core. Two cases for AE and CE in the two hemispheres that are mainly modulated by horizontal current fields are shown in Figure 7(e–j), forecasted based on the altimeter-, GLORYS12V1-, and drifter-derived current fields, respectively. The drifters provide the most accurate measurements of the near-surface ocean current velocities in the global ocean (Laurindo, Mariano, and Lumpkin 2017), providing us a valuable opportunity to access higher accurate trajectory forecasts as long as they are trapped inside eddies (Figure 7(e and f)). A visual comparison of the true and predicted trajectories (Figure 7(g and h)) shows that the predicted trajectories based on the drifter-derived currents match the true trajectories significantly better than those derived by altimeter and GLORYS12V1 ones, and to a considerable extent coincide with the advection of horizontal flow at the depth of its eddy core, as shown by the grey streamlines. The most noticeable deviations are found when using the surface geostrophic current derived from the altimeter data, while the deviations are significantly corrected when using the GLORYS12V1-derived ones. As further depicted in Figure 7(i and j), the drifter-derived forecasts yield systematically smaller errors than the other two from day 1 to day 10. Calculating the $\overline{GCD}$ for these two forecasted paths, the error using the altimeter-derived current is up to 3.53 times higher than that of the drifter-derived currents. However, the forecast accuracy is improved significantly when using GLORYS12V1-derived current which is only about 1.35 times larger than the drifter, indicating the validity of improved flow field estimation for enhancing eddy trajectory prediction. Nevertheless, with the continuous progress of current observation, including improvement of drifter coverage and the development of oceanic currents reconstruction techniques (e.g. Sinha and Abernathey 2021), further enhancements in eddy forecasting results through the integration of more precise current field are expected.

Figure 7. Case study of the effect of current fields on the eddy trajectory forecasts. (a-b) The true and forecasted trajectories of a westward propagating AE (a) before and (b) after considering the vertical shear of the background mean flow. (c-d) as in (a-b) but an eastward propagating CE. The streamlines represent the background PV gradient contours. (e and f) Eddy trajectory (grey lines) and corresponding matched drifter trajectory (dark blue lines); (g and h) Predicted trajectories for an AE (g) and a CE (h). The black and red lines are the historically modeled and the true eddy trajectories, respectively. The purple, green, and pink lines represent the eddy trajectories forecasted based on the drifter-, the GLORYS12V1-, and the altimeter-derived currents, respectively. The streamlines represent the horizontal subsurface current fields averaged over the predicted 10-day period. (i and j) The corresponding great cycle distance between the true and predicted trajectory of (g and h).

Next, we implement the medium-range prediction for over two million individual eddy tracks from 1 January 2010 to 31 May 2020 to verify the statistical reliability of the proposed forecast algorithm. The ten-day forecasted trajectories are grouped together relative to their predicted beginning positions and are presented in data density maps (Figure 8(a and b)). A striking full-angle consistency can be observed by comparing the color with the overlaid contours which correspond to the true eddy trajectories, suggesting that the performance of the proposed PTEM forecast algorithm is valid and reliable for eddies moving toward any direction. It becomes more evident by the nearly overlapping averaged predicted and true trajectories (light grey and grey lines in Figure 8(a and b)). Meanwhile, slight equatorward and poleward propagation trends can be found in both true and forecasted AE and CE trajectories, conforming to the polarity-based meridional eddy propagation, i.e. the small meridional offsets of cyclones and anticyclones are polar-oriented and equatorial-oriented, respectively (Chelton, Schlax, and Samelson 2011), exhibiting the polarity sensitivity of the proposed algorithm. Furthermore, a comparison of the hemispherical direction probabilities of eddy propagation for AE and CE in the radar plots suggests that the predicted eddy propagation can be well matched to reality (Figure 8(c and d)). Notably, the proposed algorithm’s robustness can be further verified by comparing the curves’ deviation in each direction with respect to the center. Specifically, AE has a significant southerly predominance in the northern hemisphere and a significant northerly predominance in the southern hemisphere, showing a more equatorward deflection. The eastward motion with distinct southern hemisphere preference is also evident, due to the large portion of eastward eddies in the ACC driven by the strong current system and wind stress.

Figure 8. (a–b) Density maps of the medium-range predicted eddy trajectory relative to their initial location. Overlaid contours are the corresponding data density of true eddy trajectories. The grey (light grey) lines represent the averaged true (forecasted) eddy trajectories. (c–d) Cumulative directional probability map of global eddy propagation for the predicted and true eddy trajectories in two hemispheres. The top and bottom rows are for AEs and CEs, respectively.

5. Global verifications of successively predicted eddy trajectories

Since the forecast performance of the medium-range eddy propagation has been individually and statistically proven, it is worth validating the complete trajectories of eddies over their entire lifetimes, so as to clearly display a successive, oceanwide pattern of predicted trajectories.

As an example, Figure 9 shows two cases of the successively predicted long-lived AE and CE (LAE and LCE) trajectories at 3-day, 5-day, and 10-day stages. Going through the eight panels, two characteristics become visible. First, the predicted trajectory with a 3-day stage (magenta and blue lines) is highly consistent with the true trajectories (grey lines) of LAE and LCE. As the day stage increases, some burr signals appear in the predicted trajectories, particularly at the eddy turning. The smoothness of the trajectory is an important prediction accuracy measure from a visualization perspective. The smoother the predicted trajectory, the smaller the deviation of the predicted motion from the true one, by which the forecasted eddy propagation always follows the correct trajectory. It is estimated that the averaged distance errors for the three types of stages are 5.6, 7.3, 11.4, and 5.2, 6.7, 9.4 km for the LAE and the LCE, respectively, indicating that short-range eddy motions can be predicted precisely based on the proposed algorithm, while their motions for the medium-range can be predicted accurately despite the relatively larger error. Second, the trajectory prediction is comparably accurate and stable during the maturity stage of the eddy, while tends to become relatively inaccurate during birth and decay. This is evident in the prediction error variations over the eddy’s life cycle (Figure 9(d and h)). Previous studies have noted that the submesoscale ageostrophic energy level is relatively high in eddies’ formation and decaying phases, and is low in their mature phase. Vulnerability to external deformation and larger energy variations in their formation/decaying phases (Zhang and Qiu 2018), provides a theoretical explanation for the relative poor predictability of unstable, weakly geostrophic eddy trajectories, as further discussed in Figure 13.

Figure 9. Case study of the successively predicted eddy trajectories in 3-day, 5-day, and 10-day steps. (a–c) An AE trajectory that survived 688 days from 16 June 2016 to 5 May 2018. Colored (grey) lines are predicted (true) eddy trajectories; (d) the corresponding distance error of (a–c). (e-h) same as (a-d) but is a CE trajectory that survived 831 days from 22 August 2016 to 1 December 2018.

Figure 10 exhibits the successive trajectories predicted at ten-day steps globally and regionally. Note that tracks displayed globally are randomly sampled at a ratio of 1:2 for clear presentation. A high degree of consistency is observable between the predicted and real trajectories in the global ocean, with the averaged $\overline{GCD}$ being 9.97 km. The forecast is considered precise enough to monitor eddy activities, given that the distance error is much smaller (by an order of magnitude) than the eddy scale of O (100 km). We further estimate the prediction accuracy of the westward and eastward propagating (i.e. final positions of the eddies relative to their initial locations) eddy trajectories. The averaged $\overline{GCD}$ of westward and eastward propagating eddies shows 8.33 and 11.84 km, respectively, indicating that the westward eddies are a bit more predictable than the eastward eddies. In particular, it is estimated that the 7-day averaged GCD in the SCS (100°E∼125°E, 9°S∼36°N; reference to the regional coverage of Wang et al. 2021) is 6.63 km, which is significantly smaller than the previously reported forecast error in the SCS of more than ∼13 km (e.g. Li et al. 2019; Wang et al. 2021). Moreover, the overall pattern of the enlarged typical regions with active eddy motions in different ocean basins (Figure 10(a–d)) confirms that although eddies may move randomly oceanwide (Ni et al. 2020), the proposed forecast algorithm is applicable to eddies in all directions and velocities. A careful comparison between the true and predicted trajectories is supported by the enlarged plots, indicating that most of the predicted trajectories smoothly match the true trajectory throughout the eddy’s lifecycle for both AEs and CEs. The estimated average errors in the four regions are 11.18, 9.36, 9.33, and 10.11 km for AE and 11.34, 9.28, 9.58, and 10.17 km for CE, respectively. The forecast errors in the NWP (Figure 10(a)) are larger than in the other three ocean basins because of the relatively poor forecasting performance for the eastward propagating eddies in the Kuroshio.

Figure 10. Top panel: global distribution of the true and successively predicted eddy trajectories. The red and blue (orange and green) lines are the true (predicted) trajectories for AEs and CEs, respectively. Purple boxes mark four typical areas with active eddies in the Northwest Pacific (a), the North Atlantic (b), the Southeast Pacific (c), and the South Indian Ocean (d). Lower panels are enlargements of a-d areas.

Another way to validate the predicted trajectories is to compare the eddy propagation speed derived from the predicted tracks with the altimeter-observed ones. To fully reflect the prediction performance, eddy speeds in four directions are calculated separately and shown in Figure 11. An overall impression from these four panels is that the eddy velocities calculated from the predicted trajectories are stable in all directions without much data noise, and in good agreement with the velocities derived from the real trajectories. Notably, analysis of the SSH signals in longitude–time diagrams and the observed eddy translation speed both reveal that the transition from a waves-dominated field to an eddy-dominated one occurs between latitudes ±25°∼±30° (e.g. ± 30° in Tulloch, Marshall, and Smith 2009; 25.5°S in Oliveira and Polito 2013; 27°N and 26°S in Chen, Chen, and Cao 2022), where eddies almost travel at the phase speed of the baroclinic Rossby waves. This happens to be the latitudinal zone that corresponds to the rapid improvement of the eddy forecast accuracy (see Figure 12). Quantitative calculations demonstrate that the difference between the eddy speed obtained from the predicted and the true trajectories is 0.11 km/day for the global average, further confirming the reliability of the trajectory prediction and its applicability to eddy kinematics.

Figure 11. Global distributions of eddy propagation speed derived from the predicted trajectories. Black contours are the corresponding eddy propagation speeds derived from the true trajectories. (a–d) represent the westward, eastward, northward, and southward propagation speed components, respectively.

Figure 12. (a–c) Global distribution of the eddy trajectory forecast error. From top to bottom, rows show the absolute distance error, zonal distance error, and meridional distance error, respectively. (d–f) are the zonal averages of the error corresponding to (a–c).

To perform a comprehensive quantitative assessment of trajectory prediction accuracy at a global scale, the geographical distributions and the corresponding zonal variations of the absolute, zonal, and meridional forecast error averaged per 3°×3° cell are given in Figure 12. It is estimated that the global averaged distance error is 9.89 km. Further separating Figure 12(a) into zonal and meridional components reveals the average errors of 6.47 and 6.05 km, respectively. The errors are zonally pronounced in the tropical ocean where eddies move westward without significant meridional deflections (Figure 12(b)). In contrast, larger meridional errors are evident in regions characterized by strong western boundary currents (Figure 12(c)). It can be concluded that the forecast error is significantly geographically correlated which is mainly attributed to the following factors. First, the most notable errors occur in the tropical ocean where the altimeter-based eddy detection is less effective due to both the weakened Coriolis force and the insufficient sampling by present-day altimetric technology. It is argued that the minimum altimetrically resolvable scale in the present stage stops at ∼40–60 km (Chelton, Schlax, and Samelson 2011; Pujol et al. 2016), and Chen, Chen, and Huang (2021b) have reported that more than 40% of eddies in the tropical ocean are actually undetectable by available altimeters. Meanwhile, for those altimetrically identified eddies, an unnegligible portion is misclassified in terms of polarity which is referred to as ‘abnormal eddies’. As noted in recent literature, the abnormal eddies occur in the global ocean accounting for 20–30% (e.g. Liu, Zheng, and Li 2021; Ni et al. 2021), especially in the tropical ocean within ∼15° (Chen, Chen, and Huang 2021b). This limitation and error in eddy tracking induced by the current existing observation capability fundamentally affect the performance of eddy forecasting. Second, external forcings such as the strong current and the topography affect the accuracy of trajectory forecasts. Regions with larger errors of more than 10 km are observed in the active current systems, such as the Kuroshio, the Gulf Stream, the Agulhas Current, and the Brazilian Current, whereas errors are relatively smaller in the eastern Pacific and the eastern Atlantic. Although efforts have been made to incorporate the effects of the current field into the prediction process, the prediction of complex eddy trajectories due to transient current shear (e.g. the clockwise rotating eddy trajectories near the Brazilian Current, Chen, Chen, and Cao 2022) should see further improvement as better current measurements become available in the future. Meanwhile, the error is relatively larger in the nearshore regions. It is estimated that the average prediction error reaches ∼15 km when the seafloor depth is less than 1000 m, while for depths greater than 1000 m, the average prediction error is reduced to ∼8 km. Further analysis of the relationship between topography on prediction error can be found in Figure 13(a). Third, the intrinsic eddy properties (e.g. eddy nonlinearity, eddy geostrophic ratio, lifetime) influence its predictability. Eddy motions become more challenging to predict when they are in an unstable state, with complex dynamics dominated by ageostrophic forces. For example, in the Central North Pacific, there is an unstable ‘eddy pool’ with low eddy nonlinearity (i.e. the ratio of the eddy’s maximum rotational speed and propagation speed) and high ageostrophic energy (Chelton, Schlax, and Samelson 2011; Zhang et al. 2019), making it relatively tough to forecast such unstable eddy trajectories. Therefore, although the proposed algorithm maximized the eddy forecasting accuracy, practical difficulties exist for eddy track prediction and the entire eddy remote sensing observation from identification to track. It is expected that the next-generation interferometric altimeters (e.g. SWOT, Chinese ‘Guanlan’ mission) with submesoscale ageostrophic motions resolving capability (Chen et al. 2019b; Morrow et al. 2019) will overcome such eddies’ forecasting challenges.

Figure 13. Variations of averaged eddy trajectory forecast error against (a) seafloor depth, (b) wind speed, (c) eddy geostrophic ratio, and (d) eddy normalized lifetime.

Furthermore, we examine the correspondence between the forecast error to the external forcings (seafloor topography and wind), and the eddy properties (geostrophic ratio, and life stages) (Figure 13). As previously reported, the local topography affects the propagation of eddies, particularly when encountering slopes, islands, or seamounts during their travels (Jacob, Chassignet, and Dewar 2002). This point is partially evidenced by correlating the accuracy of eddy forecasting with the ocean depth. In shallow seas, especially within 500 m, the prediction accuracy is significantly low with an average distance error greater than 15 km (Figure 13(a)). The forecast error stabilizes when depths exceed 1000 m since the influence of topography on the eddy motion becomes insignificant when the sea depth is much greater than the eddy core depth (∼226 m on the global average, Chen, Chen, and Huang 2021b). It has been noted that wind stress plays a nonnegligible role in eddy motions through eddy-wind interactions (e.g. Eden and Dietze 2009; Wilson 2016). To be specific, eddies are systematically energized when the wind stress shear is in the same direction as their rotation, while are damped if the wind stress shear is in the opposite direction (Xu, Zhai, and Shang 2016). This variation of eddy energy will inevitably affect its dynamics and propagating conditions. By estimating the averaged trajectory prediction errors at varying wind speeds, it is found that the prediction error slightly rises as the wind speed increases (Figure 13(b)). In the meanwhile, the intensity of eddy geostrophic ratio has a significant effect on forecasting accuracy. The correlation between the eddy geostrophic ratio and the $\overline{GCD}$ shows that the highest prediction accuracy is achieved when the eddy geostrophic ratio equals 1 (i.e. the eddy scale is equal to the Rossby scale, Figure 13(c)). The forecast error is significant when the eddy is under strong ageostrophic energy, but rapidly decreases when it becomes a geostrophic eddy. Finally, we try to understand how the forecast error varies during eddy’s life cycle by normalizing the error against eddy’s lifetime (Figure 13(d)). The forecast accuracy oscillates slightly at the beginning of the eddy life cycle. It gradually increases as the eddy goes towards its decaying phase, reflecting that the eddy becomes less stable with the energy gradually dropping, introducing a higher prediction error. In particular, this variation of the eddy trajectory prediction accuracy over the life cycle can be partially mapped to the geographical distribution of the forecast errors. Prediction errors are relatively small in frequent eddy formation regions along the eastern boundaries (e.g. the North Pacific and the North Atlantic), while with the energy dissipation during eddy propagation processes, the prediction accuracy becomes relatively low in the region where eddy decay is more frequent (see Figure 4(c) in Chen and Han 2019a for global distribution of the originations and terminations for long-lived eddies). Further corrections are considered necessary in the future regarding the inclusion of the above-mentioned local oceanic parameters as well as eddy ageostrophic properties and the eddy’s lifecycle in the model.

6. Discussion

Given that the prediction of eddy’s future motions should be applicable to eddies tracked based on the SLA fields that are updated in near-real-time (NRT) and completely devoid of future information, which would provide a higher scientific and practical significance. In this section, we examine the capability of forecasting eddy trajectories on the NRT SLA product with a delay of a few hours to 1 d (Pujol et al. 2016). The NRT-derived eddy trajectories (hereafter NRT eddies) are forecasted using the proposed PETM for a medium range with ten days lead time and then compared to their corresponding true trajectories tracked from the SLA in a delayed-time (DT) mode (hereafter DT eddies). The two eddy trajectory datasets are generated and quality-controlled following the methodology developed by Mason, Pascual, and McWilliams (2014) and Pegliasco et al. (2022) ranging from 1 January 2018 to 31 May 2020, which is being updated and available on AVISO. By performing a spatial–temporal matching of NRT eddies and DT eddies, a total of 3611 AE trajectories and 4289 CE trajectories with a life cycle longer than 30 days can be obtained, and the average distance between NRT and DT eddy cores is 17.69 km.

The forecast errors from day 1 to day 10 are computed in two aspects: the distance between forecasted NRT eddy trajectories and true NRT eddy trajectories (GCD_NRT), and the distance between forecasted NRT eddy trajectories to the corresponding true DT eddy trajectories (GCD_DT). It is estimated that the ten-day averaged GCD_NRT and GCD_DT are 17.01 and 24.16 km, respectively. As shown in Figure 14(a), the GCD_NRT and GCD_DT are in a range of 7.13–26.11 km and 17.05–31.59 km on 1–10 days. More importantly, as the forecast days increase, the differences between the GCD_NRT and the GCD_DT progressively decrease (from ∼10 km on day 1 decrease to ∼5 km on day 10), indicating that the eddy trajectories predicted by the NRT product tend to be closer to the true DT eddy trajectories. As can be further found in Figure 14(b–i), several individual eddies moving westward and eastward in both hemispheres are given as examples. A careful look reveals that the predicted eddy trajectories based on the NRT product are initially close to the true NRT eddy trajectories, while appearing to gradually move closer to the true DT eddy trajectories (see Figure 14(g and e) for examples). In conclusion, even when predicting eddy movement obtained from the NRT products that contain a larger intrinsic mapping error (e.g. Desai and Haines 2003), a promising predicted trajectory can be projected and it has a tendency to approach the more accurate true eddy trajectory tracked from the DT products, which not only further demonstrates the effectiveness of the proposed forecast algorithm but also proves from the side of the advantages of the DT SLA products.

Figure 14. (a) The daily-mean NRT eddy trajectory forecast errors from 1 January 2018 to 31 May 2020. The orange and pink lines represent the forecast error compared to the true NRT eddy trajectory and the true DT eddy trajectory, respectively. (b–i) Case studies of forecasting eddy trajectories moving westward (b–e) and eastward (f–i) in the two hemispheres. The first two rows are AEs, and the last two rows are CEs. The black (grey) lines are the historical NRT (DT) eddy trajectories, the pink (orange) lines are the true NRT (DT) eddy trajectories, and the green lines are the forecasted eddy trajectories based on the NRT product.

7. Concluding remarks

In analogy to atmospheric cyclones, the precise forecast of oceanic eddy trajectory is important for monitoring ocean ‘weather’ variability. Using the merged sea level anomaly products from a combination of altimeters from 2010 to 2020, an eddy trajectory forecast algorithm has been developed, which can be used precisely for medium-range eddy track prediction with a ten-day lead time. The results demonstrate the proposed methodology’s effectiveness, reliability, and robustness.

The physics-embedded temporal evolution model is constructed based on the altimeter-derived eddy propagation characteristics and the practical theory of physical oceanography. On the one hand, the integration of theoretical eddy phase speed, background current field, and eddy properties introduces novel constraints for eddy forecasting for the first time. The inclusion of these guiding factors provides an in-depth combination of the linear theory and background current fields on the modulation of the eddy propagating trajectory, resulting in a 21% improvement in forecasting accuracy $R_{Adj}^2$ and a 7 km reduction in great circle distance errors. On the other hand, the flexible modeling window and temporally evolving forecast for individual trajectories are critical strategies for realizing precise forecast, which is a significant solution to avoid errors due to manually selected historical modeled tracks, enabling more flexibility and generalizability of the proposed algorithm to permit trajectory forecasting at the global scale.

Millions of eddy trajectories have been forecasted for a medium range over the last decade, and the forecasts are highly reliable in individual cases and statistical comparisons. The result shows that the great circle distance between the eddy cores of the forecasted and true trajectories is a striking 9.9 km (8.3 and 11.8 km for westward and eastward propagating eddies, respectively) on the ten-day average, which is significant progress compared with the deep learning forecast models reported recently. By predicting the complete trajectories of eddies over their entire lifetimes, the successive patterns of predicted trajectories are presented and validated in the global ocean. We reveal that the trajectory forecast accuracy is significantly geographically correlated. Larger errors can be found in the tropical ocean of more than 12 km, while it is less than 8 km in the extratropical ocean. In particular, in areas with complex eddy trajectories and ageostrophic dominance, relatively large errors occur as well. The seafloor topography and local wind speed are also found to affect prediction accuracy. A more refined and long-range eddy trajectory forecast is expected when further factors are considered (e.g. wind stress, topography, eddy-eddy interactions) and better measurements become available in the future.

Further validation of the proposed eddy trajectory forecasting scheme on the NRT SLA product yields relatively high prediction accuracy. Impressively, the NRT-based forecast trajectories strikingly show a trend toward the real trajectories obtained from the DT product. It is expected to be integrated with eddy identification and tracking algorithms to form a complete closed-loop eddy remote sensing observation system. Moreover, it would hopefully provide a critical enhancement of various sea–air models by injecting future eddy conditions, and is of great benefit to the guidance of the marine fisheries, underwater navigation, and buoy observations.

Author contributions

Xiaoyan Chen performed methodology, software, validation, data curation, visualization, and writing of the original draft. Ge Chen contributed to conceptualization, investigation, writing review and editing, supervision and funding acquisition. Linyao Ge performed the model framework construction and validation. Chuanchuan Cao performed the investigation and validation. Baoxiang Huang performed the methodology and funding acquisition.

Acknowledgments

This research was jointly supported by the Laoshan Laboratory (No. LSKJ202201406), the National Natural Science Foundation of China (Grant Nos. 42030406 and 42276203) and the Laoshan Laboratory (No. LSKJ202204301). We thank the two anonymous reviewers for their valuable comments. We also thank Dr. Amores for providing us the helpful model eddy dataset.

Data availability statement

All data used in the analysis are available in public repositories. The all-satellite merged sea level anomaly gridded products can be downloaded from https://resources.marine.copernicus.eu/?option = com_csw%26view = details%26product_id = SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047. The eddy identification and tracking dataset can be obtained from https://data.casearth.cn/en/sdo/detail/5fa668ad1f4600005e005ba3. The near-real-time eddy dataset and the corresponding delayed-time eddy dataset are respectively the Mesoscale Eddy Trajectories Atlas Product Meta3.2exp NRT and the Mesoscale Eddy Trajectory Atlas Product Meta3.2 DT provided by AVISO at https://www.aviso.altimetry.fr/en/data/products/value-added-products/global-mesoscale-eddy-trajectory-product.html. The reanalysis current product from the GLORYS12V1 is available at https://resources.marine.copernicus.eu/. The WOA-2018 temperature and salinity gridded dataset are available at https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/, and the wind and topography data can be obtained from http://www.remss.com/measurements/wind/ and https://www.ngdc.noaa.gov/mgg/global/, respectively.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was jointly supported by the Laoshan Laboratory (No. LSKJ202201406), the National Natural Science Foundation of China (Grant Nos. 42030406 and 42276203) and the Laoshan Laboratory (No. LSKJ202204301).

Appendix A.

Further testing of the eddy trajectory forecast algorithm based on different SLA fields

A.1. Forecasting eddy trajectories derived from the numerical model SLA fields

In this Appendix, we would like to examine the eddy track forecast performance on eddies identified and tracked from the model SLA fields performed with NEMO-NATL60 model configuration during 1 July 2012–31 August 2013 in the North Atlantic Ocean. The model eddy dataset was generated by Amores et al. (2018), and we predict eddy tracks with lifetimes longer than 90 days (65 AE trajectories and 119 CE trajectories in total). As shown in Figure A1(a), the predicted trajectories are in good agreement with the true trajectories for both AE and CE. By calculating the distance error, it is found that the accuracy in predicting model eddies trajectories is even smaller than that obtained from the observations (∼8.7 km). Furthermore, the predicted trajectories moving toward different directions (i.e. westward-northward, westward-southward, and eastward) for AE and CE are shown in Figure A1(b–g), respectively. An overall impression is that the predicted trajectories of these eddies coincide well with the true trajectories, even if a certain degree of path deflection occurs during the prediction process (Figure A1(c and e)). A relatively larger prediction error is found in the eastward propagating eddies, where the predicted trajectory exhibits some degree of meridional deviation, with the forecast error of ∼9 km (Figure A1(d and g)). Nevertheless, it can be concluded confidently that the proposed eddy track forecast scheme can be applied to model simulations, which would also hopefully provide a critical enhancement of ocean models by injecting future eddy conditions.

Figure A1. Eddy trajectory forecast derived from the model eddy dataset in the North Atlantic for (a) the geographical distribution of the true and successively predicted eddy trajectories, and (b-g) three examples of the predicted trajectories propagating toward different directions for AE (b-d) and CE (e-g), respectively. Detailed information can be found in the legend.

A.2. Forecasting eddy trajectories derived from the near-real-time SLA fields

Two successively predicted AE and CE trajectories at 10-day steps during their life cycle, as well as all the predicted trajectories in the global ocean, are presented in Figure A2. It is observed that the NRT-based predicted trajectories are in good agreement with the true trajectories for both AE and CE in the global ocean, demonstrating that the proposed algorithm is a reliable method for real-time prediction of eddy trajectories.

Figure A2. (a–d) Case study of the successively predicted trajectories in 10-day steps for two AEs (a–b) and CEs (c–d). (e) Global distribution of the true and successively predicted eddy trajectories. The grey lines are the true NRT eddy trajectories, and the red and yellow (blue and green) lines are the true DT eddy trajectories and the NRT-based predicted eddy trajectories for AEs (CEs), respectively.