The numerical simulation output data are obtained from the Hybrid Coordinate Ocean Model (HYCOM; hycom.org; Bleck, 2002; Halliwell, 2004). The model domain covers the north Atlantic, the Gulf of Mexico, and the Caribbean Sea from 98° W to 20° W and from 5° S to 45° N ( Figure 1A ), at 1/25° horizontal resolution and 35 hybrid vertical layers (ATL-HYCOM 1/25). Depending on the topographic characteristics and the vertical structure, the hybrid vertical coordinates are optimized to best represent the ocean dynamics of the area. The distribution of the vertical coordinates is isobaric in the mixed layer and in very shallow areas, terrain-following at the continental shelf, and isopycnal in the stratified ocean interior (Bleck, 2002). The current configuration (ATL-HYCOM 1/25) uses initial and boundary conditions from the Global HYCOM (GLB-HYCOM), at 1/12° horizontal grid resolution and 32 hybrid vertical layers (Androulidakis et al., 2016). The atmospheric forcing is obtained from the Navy Operational Global Atmospheric Prediction System (NOGAPS) and the Navy Global Environmental Model (NAVGEM). The model performance was evaluated in Androulidakis et al. (2016) and Kourafalou et al. (2016) who showed good agreement between the model output and in situ and satellite observations. ATL-HYCOM 1/25 has also served as a Nature Run for Observing System Simulation Experiments OSSEs, designed to improve model skill, and was further evaluated in Halliwell et al. (2017). The configuration is free-running (no data assimilation) and appropriate for the process-oriented purposes of this study. Furthermore, since the model domain covers both the Caribbean Sea and the GoM, it is suitable for analyzing processes that take place in both regions as well as inter-basin exchanges (Androulidakis et al., 2021). For our analysis, we focus on a 7-year period (2010-2016) to investigate the relationship between the mesoscale eddy variability in the northwestern Caribbean Sea and the LC evolution.

In our Lagrangian analysis, we aim to detect eddy structures in the Caribbean Sea that pass through the Yucatan Channel, interacting with the predominant mesoscale features in the GoM. We employ an eddy-detecting methodology, the Lagrangian-Averaged Vorticity Deviation (LAVD), developed by Haller et al. (2016), who defined Lagrangian rotationally coherent vortices from relative vorticity. As vortices are often considered regions of vorticity maxima or minima (Haller, 2005; Haller et al., 2016), efforts have been made toward defining a vortex based on that property. Haller (2005) reviewed such definitions, according to which vortices are characterized by closed streamlines that remain invariant under Galilean transformations, such as the Q-criterion (Hunt et al., 1988), the Okubo-Weiss criterion (Okubo, 1970; Weiss, 1991), and the λ₂-criterion (Jeong and Hussain, 1995). However, these criteria do not remain invariant under rotational transformations, hence the notion of objectivity, which is necessary for a generalized vortex definition, is not satisfied. Objectivity is achieved when vortices remain the same under transformations to different frames of reference (Haller, 2005). Haller et al. (2016) overcame the difficulties of defining coherent vortices from vorticity, mainly the lack of objectivity, and showed that even though vorticity is not objective, meaning that the vorticity tensor is not invariant under changes of observers, the LAVD is, in fact, an objective field. The LAVD is defined as:

where ω is the relative vorticity, ω ¯ the spatial mean relative vorticity, and F t 0 t ( x 0 ) is the flow map, which is defined as the mapping from the initial position x ₀ at time t ₀ to some later position x(t;x ₀,t ₀) at time t= t ₀+T . According to the LAVD-based definition, the centers of rotationally coherent vortices are characterized by local LAVD maxima and the boundaries are the outermost closed convex level surfaces of LAVD. Each LAVD level set is composed of fluid particles that complete the same total material rotation relative to the mean material rotation of the whole fluid mass (Haller et al., 2016). This definition along with the objective nature of their detection method makes them robust indicators of the mesoscale eddy field in oceanic flows (Haller, 2005). Unlike geodesically detected eddies, eddies detected using the LAVD approach are allowed to experience filamentation. However, filamentation is observed to typically develop tangentially, without global breakaway. As it is based on vorticity, the method described by Haller et al. (2016) may be regarded as oceanographically more appealing than the geodesic method, which produces similar results with some caveats, such as the sensitivity to parameters during computation (Andrade-Canto et al., 2020). Tarshish et al. (2018) expanded on this issue and focused on developing an index that classifies coherent structures detected using the LAVD method into three categories based on their coherency; strictly, moderately, and weakly coherent. Using different values for convexity deficiency, they calculated the mean and median coherency indices for the coherent structures detected based on a certain convexity deficiency. The large variation of their results supported their argument on the highly sensitive nature of the LAVD method regarding the tuning of the computational parameters.

However, despite the limitations, the LAVD method is considered a robust method to identify Lagrangian coherent eddies; the study by Beron-Vera et al. (2013), with the use of concrete examples, expanded on the advantages of using an observer-independent method for eddy detection as opposed to the classic SSH-based eddy detection. Moreover, Liu et al. (2019), who quantified the material coherence of Eulerian SSH eddies and rotationally coherent Lagrangian vortices detected using the LAVD method, showed that SSH eddies are highly leaky and do not maintain their material coherence over the detection time interval, as opposed to the Lagrangian coherent vortices that were detected based on the LAVD method. They also showed that the SSH eddies may be coherent in their core, but they also may not, making SSH eddies almost indiscernible from any other body of water in the ocean, in terms of material coherence.

We applied the LAVD method in the northwestern Caribbean Sea (16° N - 22° N, 90° W - 75° W, Figure 1B ) to assess eddy coherence during the study period, using the velocity fields from ATL-HYCOM 1/25, in each of the 35 layers of the model to allow for eddy detection at depth. In our calculations, the convexity deficiency was set at 10^-2. Lower convexity deficiency led to very few coherent eddies in the region of interest. According to the study by Tarshish et al. (2018), the value of the convexity deficiency chosen in the present study (10^-2) relates to coherency indices that indicate the identification of coherent structures that are not highly leaky. Although the authors developed their methodology based on outputs from a different numerical model, they argue that it should be applicable to data from other numerical models, as well. The LAVD is calculated every 15 days and integration time T was set to 15 days, a time period that is sufficient for detecting coherent mesoscale eddies using either observational or realistic model data (Beron-Vera et al,. 2018; Androulidakis et al., 2021).

After coherence is established, the boundaries are advected to examine their pathways and physical connectivity with the GoM. The boundaries are advected as contour points using ODE45 in MATLAB and Dritschel’s curvature interpolation (Dritschel, 1988). We will be using the term “boundary” for individual coherent contours and the term “eddy” when we refer to a group of boundaries in different layers that correspond to the same structure. Eddy coherence, using the LAVD method, has only been applied in fields from realistically forced, high-resolution model in Androulidakis et al. (2021), to examine the coherence and advection of anticyclonic Caribbean eddies (Caribbean anticyclones - CarAs) at the surface, as they approached the Yucatan Channel. The calculation of LAVD at depth using horizontal velocity fields at different layers from such a model to investigate eddy coherence throughout the entire water column is a novel approach, which to our knowledge, has not been applied in previous studies.

To better track the evolution of coherent structures, 2D Lagrangian particle experiments were also performed, with releases in two different layers either inside or in the vicinity of coherent boundaries. The particles were advected using the boundary advection method that is described in the previous paragraph.

We calculate potential vorticity (PV) fluxes at the Yucatan Channel attributed to coherent structures in the Caribbean Sea as an indicator of physical connectivity between the two basins. The potential vorticity is defined as:

where f is the Coriolis parameter, z is the unit vector of the z-coordinate, u is the velocity vector from which we take into consideration only the zonal and meridional components, g is the gravitational acceleration, ρ is the potential density, and ρ0 is the reference density. The term PV ₁ represents the absolute vorticity, and PV ₂ represents the gradient of buoyancy. We compute the transport of PV (Φ) by rotationally coherent eddies as:

where P V γ ¯ is the mean PV inside the coherent boundary, A is the area enclosed by the boundary, and T is the coherence time scale. This formula was also used by Androulidakis et al. (2021) to calculate surface relative vorticity fluxes of geodesic eddy boundaries (Haller and Beron-Ver, 2013; Haller and Beron-Vera, 2014). We will discuss the PV fluxes of all the boundaries that traverse through the Yucatan Channel and into the GoM. For our analysis we will only show the fluxes attributed to PV induced by horizontal shear, which is given by: P V h s h e a r = − g ρ 0 ( ∂ v ∂ x − ∂ u ∂ y ) ∂ ρ   ∂ z , a term derived from Eq. (2).

The PV induced by vertical shear is two orders of magnitude smaller than the horizontal shear induced PV and the planetary vorticity is unrelated to changes in the LC, as discussed by Candela et al. (2002), who calculated Eulerian PV fluxes based on velocity measurements at the Yucatan Channel.

We examine eddy coherence in the northwestern Caribbean Sea with respect to LC evolution by only accounting for the coherent eddies that traverse through the Yucatan Channel ( Figure 1B ) upon 15 days of advection. Figure 2 shows the area of only the largest boundaries (diameter>=90km), normalized by the area of the largest boundary (A_max), that were detected south of the Yucatan Channel ( Figure 1B ) throughout the 7-year period along with the LC northern extension. Boundaries that are smaller than 90 km in diameter have been omitted for clarity. The dots that are assigned on the same dates usually represent parts of one eddy with coherent signature in multiple model layers, because an eddy with a coherent signature in various model layers consists of a group of coherent boundaries vertically. Thus, if multiple of these coherent boundaries have a diameter>=90km, they will represent a different dot on the same date. 86% of the coherent boundaries with diameter>=90km that are identified in the northwestern Caribbean Sea are anticyclonic, with the largest one reaching the diameter of 165 km. This boundary is part of an anticyclonic eddy that was identified as coherent in multiple layers of the model and was evident a few days before a LCE went through several detachments and reattachments prior to its final separation on August 2, 2013. Such pattern is evident in 50% of the detachments, in which anticyclonic coherent eddies with diameter>=90km were detected before LCE detachments or final separations (such as June 2010 and May 2016). In 70% of the instances where the LC was retracted, anticyclonic coherent boundaries with diameter>=90km are also evident (such as December 2011 and June 2015), while the cyclonic ones do not seem to follow a particular pattern.

Eddy coherence was evaluated throughout the entire water column (list of approximate layer depths in Table 1 ), and coherent boundaries were detected from layer 1 to layer 20 (≃650m) ( Figure 3A ). No coherent boundaries were detected below layer 20 in the time period of 2010-2016. More specifically, more coherent anticyclonic boundaries than cyclonic ones were detected in all layers (up to layer 19) except for layer 20 where only a few cyclonic boundaries were detected ( Figure 3A ). The larger number of the coherent anticyclonic boundaries in all layers suggests that coherent anticyclonic eddies in the region are more likely to have a deeper signature than the cyclonic ones. However, overall, the polarity does not play a role on the depth to which the eddies can be detected, since both cyclonic and anticyclonic boundaries are present down to ≃650m. Taking a closer look only at the boundaries with diameter >= 90km, a total of 32 large anticyclones were detected (down to ~200m), much higher than the total number of 4 cyclonic boundaries (down to ~100m), throughout the 7-year period ( Figure 3B ). Such finding indicates that large coherent anticyclones are deeper and dominate the region in comparison with the cyclonic structures of equivalent size. The actual number of 3D eddies (grouping of boundaries at different layers with similar coordinates that represent parts of an individual eddy) is not calculated and would be outside of the scope of the study, since we do not focus on the number of eddies or boundaries but on the effect of the total volume and vorticity fluxes.

Figures 4A, B show the total anticyclonic volume in each of the 15-day period (T = 15 days: LAVD assessment period) of all the structures that were detected crossing the Yucatan Channel, down to layer 20 (≃650m) and at surface, along with the LC extension, respectively. The discrete values of the anticyclonic volume of coherent boundaries are calculated at time t ₀ in the entire water column ( Figure 4A ) and at surface ( Figure 4B ). The volume of each boundary is computed by multiplying the area of the boundary by the layer thickness. The total anticyclonic volume attributed to coherent boundaries at the surface is two orders of magnitude smaller than the respective volume of coherent boundaries down to layer 20, with the pattern being similar in both cases ( Figures 4A, B ). However, coherence assessment in the entire water column provides a more comprehensive picture on the presence of LAVD vortices and their interaction with the LC system, since detection in the surface layer accounts for about 50% of the total flux occurrences. Thus, coherent boundaries that represent coherent eddies may appear at depths below the surface and would be unaccounted for if only the surface circulation was analyzed. However, it should be noted that the presence of noise in the surface layer can contribute to a failure in the LAVD method to identify vortices, hence the lower number of coherent boundaries detected in layer 1.

Focusing on the at-depth coherence assessment, the largest values of anticyclonic volume (>=10⁶ m³) are mostly evident when the LC is retracted or before LCE detachments ( Figure 4A ), indicating that large coherent anticyclones are associated with retracted LC phases (Androulidakis et al., 2021). This finding is in agreement with the result presented in Figure 2 , in which the normalized areas of the largest anticyclonic coherent boundaries are also associated to retracted LC phases. After grouping all the identified coherent boundaries based on their polarity, the anticyclonic eddies that cross the Yucatan Channel tend to also be larger than the cyclonic ones in terms of total initial volume (not shown). The 7-year mean coherent anticyclonic volume that originates in the northwestern Caribbean Sea and ends up being advected into the GoM is ~ 20 km³, whereas the mean cyclonic one is ~13 km³.

Figure 5 shows the discretized horizontal shear induced PV flux, calculated for every coherent boundary and summed over the total number of detected boundaries in all layers of the model that propagate into the GoM, along with the LC extension from 2010 to 2016. The PV flux is summed over the 20 layers of the model in which coherence was detected and it is computed within the boundaries that cross the red line in Figure 1B . The coherent PV flux values are normalized by the maximum value of the examination period. The largest coherent PV fluxes are anticyclonic (negative values) and, in 5 cases, intensified anticyclonic PV fluxes are evident before LCE detachments or final separations (such as June 2010 and December 2012, Figure 5 ). This suggests that the presence of anticyclonic eddies in the northwestern Caribbean Sea is connected to the LC evolution and more specifically to the LCE detachments or separations. Furthermore, large anticyclonic fluxes are evident during retracted LC phases ( Figure 5 ) and along with the findings from Figure 4 indicate that the coherent anticyclonic activity in the northwestern Caribbean Sea promotes retracted LC phases, in agreement with Androulidakis et al. (2021). The Eulerian study by Oey et al. (2003) showed that anticyclonic activity in the Caribbean Sea is also associated with retracted LC phases, although a direct comparison to our results would not be optimal due to the lagrangian nature of our study. On the other hand, Candela et al. (2002), in their Eulerian study, showed that cyclonic vorticity flux is present before LCE detachments. In our study, there are a couple of cases (i.e., 2015, 2016) out of a total of 9, where cyclonic activity is pronounced right before the start of LCE shedding events ( Figure 5 ). Nonetheless, the values are not as large as the largest anticyclonic vorticity flux values. However, larger cyclonic PV fluxes (>0.2) do occur in other instances (years 2013 and 2014) as well, when the LC is extended, about 45-50 days before LCE detachments/separations. This is a considerably larger time interval than the coherence assessment time interval selected in the present study (15 days), thus establishing a relationship between intensified cyclonic PV fluxes and LCE detachments is not trivial. The controversy in the results as well as the high variability of the system indicated that a statistical approach using a longer time series of data may be needed to further elucidate the relationship between the coherent eddy activity in the Caribbean and the LC evolution.

During extended LC phases, ~35% of the normalized coherent PV flux values are in the bin range between 0 and 0.1, indicating weak coherent mesoscale activity in the northwestern Caribbean ( Figure 6 ). The rest of the occurrence frequencies during extended LC phases are equally distributed in positive and negative PV flux values, except for a 7% of increased negative PV flux values occurring only when the LC is extended (north of 26° N). This is possibly connected to intensified anticyclonic activity before LCE detachments, as presented in Figure 5 . On the other hand, during retracted LC phases (south of 26° N), almost 70% of the coherent PV flux values are negative, a finding that agrees with Androulidakis et al. (2021), who argued that anticyclonic surface relative vorticity attributed to geodesic eddies was more pronounced when the LC is retracted.

We now focus on specific case studies of coherent Caribbean anticyclones associated with different LC phases.