Mechanisms of the overturning circulation in the northern Red Sea other than convective mixing

The northern Red Sea (NRS) is where Red Sea Outflow Water (RSOW) and, occasionally, Red Sea Deep Water are formed. Glider observations are used to describe the formation mechanisms and pathways of intermediate waters in the NRS in late winter from 31 January to 18 April 2019. Utilizing glider observations, atmospheric reanalysis products, and satellite datasets, we evaluated the mesoscale activity and atmospheric conditions that contribute to outflow water formation. The cyclonic circulation in the regional surface dense water exposes it to the atmosphere, ventilating the water column and contributing to phytoplankton growth – enhancement of chlorophyll concentration due to nutrients upwelled into the euphotic layer (Z_eu). Subduction of this water in the three-dimensional cyclonic circulation transported oxygenated and elevated chlorophyll concentration water to depths between 150 and 250 m along the 28.2 kg m⁻³ isopycnal. Unlike previous observations, in late February, a strong anticyclonic circulation blocked the inflow of warmer, fresher water into the region. It was accompanied by a negative heat flux and uplifting of dense water to the surface. Net cooling through mid-March cooled the incoming surface waters from the south. At the end of the observational period, the intrusion of warmer, fresher waters from the south coincided with the re-establishment of cyclonic circulation and capped the dense surface water that had formed during March. These observations demonstrate that multiple processes contribute to RSOW formation: convective mixing, cyclonic uplifting of dense water, subduction, and mesoscale or submesoscale processes.

1 Introduction

The Red Sea (RS) is a narrow, elongated, and meridionally oriented basin lying between the Asian and African continents. Its subtropical location results in significant buoyancy losses due to high evaporation (nearly 2 m yr⁻¹), negligible precipitation, and effectively no riverine inputs, resulting in it being one of the world's saltiest and warmest seas (Edwards and Head, 1987; Smeed, 1997, 2004; Sofianos et al., 2002). The RS experiences seasonally reversing winds over its southern region, coupled with the Arabian Sea's monsoonal forcing. The reversing winds control the circulation and water mass exchange through the Strait of Bab al Mandab (Abualnaja et al., 2015; Patzert, 1974). These processes, along with the buoyancy forcing, drive the large-scale circulation (Bower and Farrar, 2015; Murray and Johns, 1997; Patzert, 1974; Sofianos and Johns, 2007). The inflow of water from the Gulf of Aden compensates for the evaporative water loss in the RS. The advected northward waters contribute to the overall latitudinal gradient in salinity and temperature from south to north. The northward advection of comparatively fresh and warm water from the south has an important role in the stratification (Asfahani et al., 2020; Churchill et al., 2014; Sofianos and Johns, 2007; Zarokanellos et al., 2017b), with the annual flux into the RS reaching up to 0.22 Sv (Sofianos and Johns, 2002). This water has been traced to the NRS at 28° N (Zarokanellos et al., 2017a). It is modified through heating, evaporation, and mixing as it progresses northward (Cember, 1988; Sofianos and Johns, 2003; Sofianos and Johns, 2007).

Significant mesoscale activity is found along the RS's main axis (Morcos, 1970; Morcos and Soliman, 1972; Quadfasel and Baudner, 1993; Sofianos and Johns, 2007; Zhan et al., 2014), which results from baroclinic instabilities (Zhan et al., 2016) and the presence of the Eastern Boundary Current (EBC; Zarokanellos et al., 2017a, b; Biton et al., 2008, 2010; Bower and Farrar, 2015; Eshel et al., 1994; Eshel and Naik, 1997; Sofianos and Johns, 2007). Some of these mesoscale eddies are thought to be quasi-permanent features, in particular the cyclonic eddy (CE) in the northern Red Sea (NRS) and the anticyclonic eddy (AE) in the central Red Sea (CRS) (Chen et al., 2014; Yao et al., 2014b). The CE in the NRS plays a crucial role in forming the RSOW (Asfahani et al., 2020; Sofianos and Johns, 2007). Zarokanellos et al. (2017a) showed that an AE in the CRS can redirect or deflect the advected northward flow of water from the Gulf of Aden. These mesoscale eddies substantially affect heat and salt advection and the distribution of biogeochemical properties (Chen et al., 2014; Raitsos et al., 2013; Triantafyllou et al., 2014), which are often useful tracers of physical processes. Eddies also transfer energy and momentum to the mean flow, driving the circulation (Lozier, 1997). In addition, mesoscale features fundamentally modulate local nutrient fluxes and phytoplankton dynamics (Longhurst, 2007). Eddy activity in the RS has shown that it significantly influences optical and biological properties (Kürten et al., 2016; Mahadevan, 2015; McGillicuddy et al., 1998; Pearman et al., 2017). In particular, mesoscale eddies in the RS facilitate the vertical transport of nutrients into the upper water column, thereby impacting biogeochemical variability (Zarokanellos et al., 2017b; Kürten et al., 2019).

Mesoscale processes can induce subduction events over larger areas with significant impacts on the marine ecosystem (Mahadavan et al., 2016). Subduction is the process where physical and biogeochemical tracers are transported from the surface to the interior. This, in turn, influences the chlorophyll (CHL) flux in Z_eu, playing a significant role in the RS's biological variability. Furthermore, while the upper ocean is primarily ventilated through subduction events, the ventilation in the deep ocean occurs primarily through open-ocean convection (Asfahani et al., 2020; Williams and Meijers, 2019). Seasonal mixing, eddy interactions, and cross-shelf exchange drive biogeochemical fluxes in the upper layer. In the RS, seasonal phytoplankton variability has been studied primarily through remote sensing (Acker et al., 2008; Gittings et al., 2018; Racault et al., 2015; Raitsos et al., 2013), with only a few studies incorporating in situ observations (Kheireddine et al., 2020; Zarokanellos et al., 2017a, b; Zarokanellos and Jones, 2021). These studies highlight ecological connections between coastal and open-sea regions, which are vital for coral communities (Acker et al., 2008; Raitsos et al., 2017). However, remote sensing is limited to surface observations, as ocean color data do not capture subsurface phytoplankton distributions. Integrating bio-optical data from both remote sensing and in situ measurements has provided valuable insights into the RS's biogeochemical dynamics (Gittings et al., 2018; Kheireddine et al., 2020; Racault et al., 2015; Raitsos et al., 2013; Zarokanellos et al., 2017a, b; Zarokanellos and Jones, 2021).

Water mass transformation typically occurs in the surface ocean in specific regions where air–sea interaction acts powerfully on the upper layer (Emery, 2001; Iselin, 1939). Significant surface cooling results in buoyancy loss that can contribute to deep mixing and convection. The NRS has been considered the main area of RSOW formation and, occasionally, Red Sea Deep Water (RSDW) formation (Papadopoulos et al., 2015; Sofianos and Johns, 2003, 2007; Yao et al., 2014b), due to the high evaporation rates and significant surface water cooling that occur during winter (Papadopoulos et al., 2013). These preconditions, in addition to the presence of the CE, where shallowing of the isopycnals occurs at the eddy center, contribute to the formation of the aforementioned RSOW and RSDW (Abualnaja et al., 2015; Asfahani et al., 2020; Sofianos and Johns, 2003; Yao et al., 2014b; Zhai et al., 2015). The newly formed RSOW gradually sinks until it reaches an equilibrium density of nearly 27.5 to 27.7 kg m⁻³ throughout the basin (Zhai et al., 2015) and flows southward, where it exits the RS through the Strait of Bab El Mandeb (Cember, 1988; Sofianos and Johns, 2003; Yao et al., 2014b). Once it enters the Gulf of Aden, it is mixed due to the intense mesoscale eddy activity at depths between 400 and 1000 m (Bower and Furey, 2010). Despite significant dilution of the thermohaline properties, RSOW has been observed as far south as the Agulhas Current below 32° S (Beal et al., 2000) and as far east as the Bay of Bengal (Jain et al., 2017).

Numerical simulations and a few in situ observations suggest three potential mechanisms for RSOW formation. The first mechanism, open-ocean convection, is associated with the presence of a cyclonic gyre in the NRS (Sofianos and Johns, 2007; Papadopoulos et al., 2015). Strong atmospheric forcing in the region and the presence of a cyclonic gyre create favorable conditions for convection events. The second mechanism is mixed layer deepening during winter, resulting from a large negative heat flux (Zhai et al., 2015). The third mechanism is the combined effect of the cyclonic gyre and the weakening of the stratification that results from strong atmospheric forcing (Chen et al., 2014; Clifford et al., 1997; Manasrah et al., 2004; Morcos and Soliman, 1972; Yao et al., 2014b). Concurrent with the formation of the RSOW, the RSDW forms in the NRS when strong cooling and evaporation occur in the gulfs of Suez and Aqaba (Cember, 1988; Papadopoulos et al., 2015; Sofianos and Johns, 2003; Yao et al., 2014b). The RSDW is evident below 300 m, while the RSOW is typically found between 200 and 300 m and can be identified by the relatively high oxygen concentration, which can be distinguished from the RSDW. The regionally formed RSOW has a significant role in the overturning circulation, as indicated by in situ observations and numerical simulations (Papadopoulos et al., 2015; Sofianos and Johns, 2003; Zhai et al., 2015). In addition, RSOW contributes to the ventilation of the RS (Papadopoulos et al., 2015; Sofianos and Johns, 2007; Woelk and Quadfasel, 1996; Yao et al., 2014b; Zhai et al., 2015) and the salt budget of the Indian Ocean due to the high-salinity concentrations present at intermediate depths (Beal et al., 2000).

The primary objective of this study is to understand the mechanisms that contribute to the mass formation of RSOW in the NRS and the biogeochemical responses associated with these physical processes. The second objective is to evaluate how atmospheric forcing affects mesoscale dynamics in the study area. This paper is organized as follows: in Sect. 2, the in situ, reanalysis, and remote sensing observations are presented. Section 3 describes the coupling of the atmospheric forcing and the in situ observations. A comprehensive description of the flow variability and the mechanism of RSOW formation is described, including the coupling between physical and biogeochemical processes in the study area. Finally, Sect. 4 presents a discussion and comparison with previous regional studies, along with the conclusions.

2 Data and methods

2.1 Glider observations

A sustained glider line that was traversed in a 3–4 d period in the NRS was used to capture the wintertime evolution of physical and biogeochemical characteristics in the NRS (Eyouni et al., 2025). The Seagliders (hereinafter “gliders”) were equipped with a conductivity–temperature–depth (CTD) sensor, a dissolved oxygen (DO) sensor, and a triplet fluorometer (Table 1). The glider was deployed along a transect oriented approximately perpendicular to the coastline offshore from Duba, referred to as the “Duba line” (Fig. 1a). The line extended from approximately 5 km off the coast to nearly 75 km offshore. The glider completed each transect in about 3.5 d and was programmed to dive from the surface to ∼750 m. The average horizontal speed was ∼25 cm s⁻¹, and each dive cycle took ∼3 h to complete, depending on the target depth, topography, and sea conditions. The deployment spanned the period from 31 January to 21 April 2019 (Table 1, Fig. 1b).

Table 1 Summary of the measured variables and their corresponding sensors deployed on the glider.

Figure 1 (a) Topographic map of the RS. The two black dashed boxes in the northern and southern RS show the study area in the north and the region in the south where exchange with the Gulf of Aden occurs. The yellow line inside the northern box indicates the glider line location, and (b) the yellow sawtooth line represents the undulating trajectory of the glider superimposed on the bathymetry, which has been reproduced from the GEBCO_2021 Grid, GEBCO Compilation Group (2021; https://doi.org/10.5285/c6612cbe-50b3-0cff-e053-6c86abc09f8f).

Each vertical profile has been examined for spikes and outliers before further analysis. No evidence of thermal lag has been found in the CT observations. The WETlabs EcoPuck Triplet Fluorometer (FL3) and Backscatter (BB3) sensor were factory-calibrated, and dark counts were measured prior to deployment in March 2019. The pre-deployment dark counts for FL3 were consistent with the factory dark counts. Roesler et al. (2017) performed a global comparison of fluorometer and extracted CHL measurements and recommended that the factory-calibrated CHL should be divided by 2. Thus, the CHL is divided here by 2 to correspond to in situ CHL concentrations as described in Roesler et al. (2017). As the RS is a region with high irradiance, the CHL can often experience quenching. In this study, we examined the vertical profiles of CHL fluorescence for quenching, and no significant quenching was detected. The DO measurements were adjusted by a correction factor based on the median oxygen saturation in the upper 10 m. The correction factor was the ratio of the median oxygen saturation concentration for the upper 10 m to the median measured concentration within the upper 10 m of the water column. All of the water column data were then corrected by multiplying the reported DO concentrations by the correction factor. The correction provided consistency between the current glider dataset and the vertical distribution for this region in the World Ocean Atlas (Zarokanellos and Jones, 2021; Garcia et al., 2019a; Johnson et al., 2015). To facilitate the data processing, plotting, and analysis, the corrected data were projected onto a grid with a vertical resolution of 2 m and a horizontal resolution of 2.5 km, which is the nominal distance between glider surfacing values during the mission. This study focuses on the upper 500 m of the water column. Mixed layer depth (MLD) was calculated for each vertical profile based on a density difference of ≤0.03 kg m⁻³ relative to the density at 10 m depth (de Boyer Montégut et al., 2004). The Brunt–Väisälä frequency (BVF) was computed to measure the stratification. Geostrophic velocity and potential density were also calculated using the TEOS-10 toolbox (TEOS et al., 2010; Lea et al., 2015). The level of no motion in the geostrophic velocity has been examined, and no significant motion below 500 m has been observed.

2.2 Remotely sensed data

2.2.1 Sea level anomaly (SLA)

Previous satellite observations and numerical simulation studies suggest a cyclonic gyre located in the NRS concurrent with convective mixing where the RSOW mass is formed (Sofianos and Johns, 2002; Zhai et al., 2015). Therefore, in this study, SLA data were used to characterize the spatial and temporal evolution of the large-scale circulation patterns of sea level and geostrophic velocity during the wintertime. The SLA data are based on the multi-mission altimeter Archiving Validation and Interpretation of Satellite Oceanographic (AVISO) data provided by the Copernicus Marine Environment Monitoring Service (CMEMS). The data have been gridded at a regular 0.258°×0.258°. The obtained daily measured data were taken within the domain from 20 to 28° N and from 32 to 42° E, inclusive of the glider deployments. Within the NRS subdomain (26–28° N, 33–37° E), the data were temporally averaged to approximately fit the periods of the glider transects. Previous work has validated the use of the AVISO product from comparisons with both in situ and numerical model results (Hernandez and Schaeffer, 2001; Zhan et al., 2014; Zarokanellos et al., 2017b).

2.2.2 Sea surface temperature (SST)

SST data were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) imagery that provides a daily image with a spatial resolution of 4 km. The data were obtained within the domain from 24 to 28° N to identify and determine the spatial and temporal evolution of the large-scale patterns of the SST (Werdell et al., 2013).

2.3 Atmospheric reanalysis products

The atmospheric parameters for this study were derived from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) reanalysis dataset (Rienecker et al., 2011; Gelaro et al., 2017). MERRA-2 offers a high-resolution representation of atmospheric conditions, providing daily mean values on a 0.5°×0.652° grid with a temporal resolution of 1 h. This dataset was chosen due to its demonstrated accuracy in representing heat fluxes within the region of interest, as validated by recent studies such as Al Senafi et al. (2019). For this analysis, data were extracted and spatially averaged over a defined box encompassing the NRS (NRS 26–28° N, 33–37° E; Fig. 1a), excluding land coverage. Daily means were calculated from the hourly data for the selected parameters wind speed and direction, surface net heat flux (Q_net), and evaporative heat flux, which is considered to be the primary driver of heat loss in the RS (Sofianos and Johns, 2002).

Q_net, produced using the Coupled Ocean-Atmosphere Response Experiment (COARE 3.0) formulation (Fairall et al., 1996), consists of the sum of shortwave (Q_sw, absorbed), longwave (Q_lw, emitted), latent (Q_L, evaporative), and sensible (Q_S, (Q_S, conductive) heat fluxes, as all of these terms are positive when they heat the water column. The quantity Q_net is calculated following formula:

$$\begin{array}{}\text{(1)}& Q_{\mathrm{net}}=Q_{\mathrm{sw}}+Q_{\mathrm{lw}}+Q_{\mathrm{S}}+Q_{\mathrm{L}}.\end{array}$$

2.4 Empirical orthogonal function (EOF) analysis

To compare the SLA from 2019 with observations from the preceding years, EOF analysis is carried out on the SLA data to evaluate the spatial and temporal patterns of the variability of the SLA data in the winter–spring period. Each eigenvector describes the spatial pattern (modes) of that variability for 5 months from January to May for the years 2016, 2017, 2018, and 2019. Only the first mode of the EOF is used to compare the seasonal evolution of the spatial pattern each year. In addition, the explained variance with the eigenvalue provides the relative contribution of a specific mode to the variability (Zhang and Moore, 2015).

2.5 One-dimensional mixed layer model

The Price–Weller–Pinkel (PWP) mixed layer model (Price et al., 1986) was applied to evaluate the local atmospheric effects on the ocean mixed layer. The model input includes the following terms: radiative heat flux (shortwave and longwave), latent and sensible heat, freshwater flux (evaporation [E] and precipitation [P]; [E−P]), and wind stress components (τ_x and τ_y). As the RS is sandwiched between two extreme desert regions, precipitation is considered to be negligible (P=0). The one-dimensional PWP model was applied to estimate the local, atmospherically driven evolution of the mixed layer depth during the observation period. The model was executed for three sequential subsets delineated from the glider observations: the cooling phase; the cool, salty, and dense phase; and the warming–freshening phase. The model was initialized with the average temperature and salinity profile for one complete glider transect at the beginning of each simulation period (30 d) and then stepped forward in 24 h (1 d) increments subject to the heat, freshwater, and momentum fluxes. The daily time step was selected as insignificant diurnal variability was observed in the mixed layer temperature and salinity from the glider.

The PWP model produces a mixed layer through a vertical exchange process between the air–sea interaction and vertical mixing. It assimilates time series of surface heat flux, wind, and precipitation and applies these forcing parameters to the initial vertical profile of temperature and salinity. The model interpolates the momentum components driven by winds, cooling, and evaporation to induce convective instability, entrainment from the pycnocline, and a mixing term generated by vertical current shear. In our case, since P=0, only surface heating affects re-stratification.

The convective adjustment in the PWP model starts with grid cells with unstable stratification being homogeneously blended with neighboring cells. The convective correction follows the bulk mixed layer parameterization, where the mixed layer deepens when the bulk Richardson number, R_ib, falls below a threshold value of 0.65 (Price et al., 1986).

The bulk Richardson number is expressed as

$$\begin{array}{}\text{(2)}& R_{\mathrm{ib}}={\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }g}{{\mathit{\rho }}_{\mathrm{0}}(\mathrm{\Delta }U{)}^{\mathrm{2}}}},\end{array}$$

where Δρ (density) and ΔU (velocity) are the differences between their values within the mixed layer and their values below the mixed layer, respectively (Price et al., 1986). The variable ρ₀ is the reference density and g is the gravitational acceleration. Then, the model adds local shear instability below the mixed layer, where mixing due to strong shear is parameterized based on a threshold gradient Richardson number, R_ig, defined as

$$\begin{array}{}\text{(3a)}& R_{\mathrm{ig}}={\displaystyle \frac{N^{\mathrm{2}}}{S^{\mathrm{2}}}}=\mathrm{0.25},\end{array}$$

where S² is the square of the shear below the mixed layer and

$$\begin{array}{}\text{(3b)}& N^{\mathrm{2}}={\displaystyle \frac{g}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{\mathrm{\Delta }z}}.\end{array}$$

The vertical resolution of the model depth bin was set to 2 m for alignment with the gridded glider data. The momentum (horizontal diffusivity) and vertical diffusivity were set to 10⁻⁵ and 0 m² s⁻¹, respectively (Sanikommu et al., 2020; Zhai et al., 2015). The maximum depth of the PWP experiment for the run and the initial depth range for the profiles of salinity and temperature were 400 m, the depth at which the divergence of hydrographic variables between summer and winter was minimal. Estimation of the MLD at each time step throughout the PWP model run used the same criteria as for the glider data by determining the depth range over which the density increase relative to 10 m was no more than 0.03 kg m⁻³ (de Boyer Montégut et al., 2004).

3 Results

3.1 Atmospheric forcing

Regional atmospheric forcing is a major factor affecting the seasonal variability of the RS circulation. The wind direction in the northern part of the RS is predominantly from north-northwest (NNW; Fig. 2a). The net heat flux was initially negative, with heat flux losses of up to 300 W m⁻² in January and 250 W m⁻² in February, transitioning in late March from net negative to net positive and thus beginning the onset of the seasonal heating period (Fig. 2b).

Figure 2 (a) Wind vectors for the NRS and (b) net heat flux for the NRS (average for the black dashed box over the NRS, Fig. 1) for the period from 1 January through 1 July 2019. The red box indicates the period of the glider mission.

3.2 Upper-ocean response to atmospheric forcing

3.2.1 Regional response from a remote sensing perspective

Remotely sensed ocean imagery demonstrates the seasonal evolution of the upper ocean during the glider deployment. Figure 3 provides images of 8 d composites for SLA and SLA-derived geostrophic velocity, sea surface temperature, and CHL concentration. At the onset of the glider deployment in late January–early February, the eastern boundary coastal flow was northward (Fig. 3a). Consistent with this northward flow, a tongue of warmer, low-CHL water extended northward along the Saudi coastline (Fig. 3b and c). These observations are consistent with previous observations from 2016 (Asfahani et al., 2020) and with the mean structure typically observed for SST in the winter months (Karnauskas and Jones, 2018). Cooler higher-CHL water is observed on the western side of the basin, perhaps due to the convective mixing described by Kheireddine et al. (2020).

Figure 3 The 8 d averages of the sea level anomaly and geostrophic velocity (a, d, g, j) from AVISO as provided by the Copernicus Marine Environment Monitoring Service (CMEMS), the 4 μ nighttime SST from the MODIS Aqua satellite (b, e, h, k), and CHL from the MODIS Ocean Color Instrument (OCI) (c, f, i, l) during the presence of the EBC (a–c), AE (d–f), pair of anticyclonic eddies (g–i), and lateral advection (j–l). The red line indicates the location of the Duba glider track. The blue box in panel (b) indicates the region of the northern RS (NRS) that was used for regional averages (∼200 km × 200 km). The dashed blue line separates the eastern half of the NRS from the western half.

Following the initial phase of northward coastal flow, an AE developed in the northeastern RS (Fig. 3d). During this period, the flow was southward across the glider line, and there was no indication of northward advection of warmer, low-CHL water (temperature > 24 °C and CHL < 0.1 mg m⁻³) from the south (Fig. 3e and f). The AE appeared to block the warmer, low-CHL water transport into the region. The surface temperature of the NRS became cooler, reaching a mean surface temperature near 22.5 °C, and the temperature difference between the western and eastern sides of the northern region decreased to less than 0.25 °C (Figs. 3h and 4a).

Figure 4 Time series of near-surface (6 m) and deep (500 m) (a) temperature, (b) salinity, (c) density, and (d) MLD for the glider deployment from 31 January through 21 April 2019. This is the complete time series, irrespective of the glider's location along the transect. In panel (a), temperature averages from the 8 d MODIS SST are shown for the entire NRS (solid blue line), the eastern half of the northern RS (NRS-East, dashed blue line), and the western NRS (NRS-West, dotted blue line). The geographical boundaries of these subregions are shown in Fig. 3b.

In the latter half of March, two anticyclonic eddies (Fig. 3g), between 22 and 26 °N, apparently blocked the northward flow of water from the south and isolated the northern part of the RS from the inflow of the warmer, low-CHL water. Consequently, the near-surface area in the NRS became almost thermally homogeneous, with small temperature variations in its spatial distribution (Fig. 3h). The densest surface water was observed in the NRS during this period (Fig. 4c). By the beginning of April, these two eddies had dissipated, and warm, low-CHL water again advected into the NRS along the eastern coastline (Fig. 3–l), re-establishing the temperature differential between the eastern and western halves of the NRS (Fig. 4a–c).

3.2.2 Upper-layer variability

As the atmosphere progressed through its typical annual cycle, the near-surface ocean demonstrated high variability. In order to compare the data accurately, the depth of 6 m was chosen to represent the near-surface layer, while the depth of 500 m was selected to represent the near-bottom layer, as it is the most isolated from the surface influence and shows the least variability within the dataset. The time series of surface (6 m) and 500 m values for temperature, salinity, and density are shown in Fig. 4. This is a continuous time series of glider data, irrespective of its location along the transect. Distinct phases are evident in the time series. The early phase, consistent with the atmospheric forcing, shows a general cooling trend from mixed layer temperatures near 24 °C in the early period to about 22.5 °C during the coolest phase, after which temperatures rose again to nearly 24 °C. For most of the period, the western half of the northern RS was cooler than the eastern half, where the glider was operating (Fig. 4a). However, during the coolest period, the temperatures were nearly homogeneous across the entire northern RS. Correspondingly, mixed layer salinities rose from values of 40–40.2 during the early phase to nearly 40.4 during the cool, salty period and then returned to values between 40 and 40.2 in the warming period. As a result, the near-surface density anomaly initially increased from 27.6–27.8 to 28.3–28.4 kg m⁻³ during the dense period. As expected, no effects of the surface forcing were apparent at 500 m.

In the latter part of the cooling phase (between 17 February and 12 March), the MLD exceeded 100 m much of the time. Although the deepest mixed layer occurred during the cool, salty phase, the MLD during this period was highly variable (Fig. 4d). The variability in the MLD is a result of the new, well-formed cyclonic eddy during that time. A shallow MLD can be present in the center of a newly formed cyclonic eddy, as is shown in Fig. 7, and it is also consistent with the patchiness of the SST (Fig. 3e).

Individual glider sections provide insight into mesoscale processes occurring in the NRS during this winter–spring period. In the early cooling phase, the isopycnals were tilted strongly downward from offshore to nearshore, with cool, dense water near the surface offshore and warmer, fresher water near the surface nearshore (Fig. 5a and b). The overall geostrophic flow was northward, with an average velocity of 0.31 m s⁻¹ in the upper 100 m (Fig. 5e), consistent with the geostrophic velocity calculated from the sea level anomaly (Fig. 3a). The maximum upper-layer stratification occurred in the 20 km offshore from the transect, where uplifting of the pycnocline resulted from the cyclonic circulation (Fig. 5d). The maximum near-surface CHL concentration occurred where dense water (≥28.1 kg m⁻³) offshore shallowed to within 50 m of the surface.

A small-scale cyclonic feature centered about 43 km offshore was embedded within the larger-scale flow (Fig. 5a, b, and e). This feature was not observed in either the previous or subsequent transects, and each was separated from the current transect by approximately 3 d. Thus, this was a transient feature on the glider line, which we conclude was advected across the glider line within the larger-scale flow. A feature characterized by elevated CHL ($\sim \mathrm{30}\times {\mathrm{10}}^{\mathrm{3}}$ mg m⁻³) and DO (∼178 µmol kg⁻¹) was observed along the 28.2 isopycnal between 20 and 40 km offshore, compared to the rest of the transect, where the CHL was $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ mg m⁻³ and the DO ∼94 µmol kg⁻¹, respectively, while a low BVF is observed in the same area. The feature spanned a depth range from approximately 150 m at 40 km offshore to 250 m at 20 km offshore, suggestive of subduction of denser near-surface water and downward transport along the isopycnal below the mixed layer and the euphotic zone (Fig. 5c, d, and f). This signal was also present in the backscatter at 650 nm (Fig. 5h). While the higher concentration of the backscatter is evident offshore at the surface ($\sim \mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹ sr⁻¹) and decreases with depth and proximity to the shore in the same area as the elevated CHL and DO, the backscatter reaches values of $\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹ sr⁻¹, whereas in the surrounding waters it is less than $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹ sr⁻¹. This aligns with subduction of CHL (Fig. 5g) and backscatter at 650 nm (Fig. 5h) along the isopycnals between 28.2 and 28.3 kg m⁻³ below the mixed layer and Z_eu (120 m) at depths greater than 200 m, contributing to the export of organic matter, as has been observed in other regions of the global ocean (Zarokanellos et al., 2022; Zarokanellos and Jones, 2021; Erickson et al., 2016).

Figure 5 Glider section of the (a) depth-averaged current (DAC), (b) geostrophic velocity, (c) potential temperature, (d) salinity, (e) DO, (f) Brunt–Väisälä frequency, (g) CHL, and (h) backscattering coefficient at 650 nm for 5–9 February 2019. Density isopycnals are shown in panels (b)–(h) with solid grey lines, while the isopycnals of 28.2 and 28.3 are distinguished by black lines (or white for panels f–h). The solid magenta line represents the MLD. Geostrophic velocity, calculated relative to 500 m, is positive north-northwest and negative south-southeast parallel to the coastline.

The glider observations of CHL, backscatter at 650, and DO allow us to independently trace the subduction with three different bio-optical tracers. Indeed, the observed elevated CHL – commonly associated with phytoplankton growth in the literature – primarily only occurs offshore, within the Z_e, where dense water (≥28.1 kg m⁻³) rose to a depth shallower than 50 m, bringing up nutrients from deeper layers. Also, the Z_eu is located around 120 m in the RS, and light at greater depths is too low to sustain photosynthesis (Zarokanellos and Jones, 2021). Furthermore, this transient eddy about 43 km offshore was not observed in either the previous or following sections, and it was embedded within the larger-scale flow (Fig. 5a, b, and e). The observed high DO concentration on the surface can be a result of photosynthesis. The co-occurrence of high CHL and DO at depths below the Z_eu suggests that this water was originally at the surface before it was transported and subducted deeper. The fact that the high-CHL and DO waters align along the 28.2 isopycnal (Figs. 5, 11a and b) indicates that their subduction is associated with an eddy wherein the denser surface water is forced below, with the lighter water following the 28.2 isopycnal rather than being mixed vertically.

Following the initial period of northward transport, the circulation changed significantly in late February, reversing the direction of flow across the glider line due to the presence of an AE in the northeastern RS (Figs. 3d–f and 6). Based on the glider sections, the southward coastal flow began in mid-February and persisted for approximately 3 weeks. The isopycnal structure associated with the anticyclonic geostrophic flow is evident below the mixed layer. In the glider section from 1 to 5 March, the geostrophic velocity varied between weakly northward at the offshore and inshore ends of the line to a maximum southward velocity of 0.67 m s⁻¹ and an average southward velocity of 0.14 m s⁻¹ in the upper 100 m, which is smaller in magnitude than the northward flow during the preceding cyclonic period (Fig. 6b). During this period, the upper boundary of the pycnocline, 28 kg m⁻³, which had been near the surface offshore during early February, was now at nearly 150 m depth in the offshore portion of the transect but rose to the surface near the coast (Fig. 6). The MLD was consistently deeper than 100 m in the offshore 30 km of the transect, and in the nearshore 40 km it varied between ∼25 and 120 m. CHL was uniformly distributed within the mixed layer with diel patterns in concentration, which appear as spatial patterns in the 3.5 d transit of the line (Fig. 6g).

Figure 6 Same as for Fig. 5 but for the period 1–5 March 2019.

Typically, the oxygenated waters are located in the surface layers within the MLD. However, the observed bolus indicated that highly oxygenated waters had been trapped below the MLD between the 28.2 and 28.3 kg m⁻³ isopycnals at 150 to 250 m depths and between 20 and 50 km offshore. The average DO concentration within the bolus is ∼177 µmol kg⁻¹, while the CHL is around $\mathrm{4.6}\times {\mathrm{10}}^{-\mathrm{3}}$ mg m⁻³. The surrounding waters below the 28.3 isopycnal indicate that the DO and CHL values reach 62 µmol kg⁻¹ and $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{3}}$ mg m⁻³, respectively. Above the 28.2 isopycnal, the DO and CHL have values of 203 µmol kg⁻¹ and $\mathrm{79}\times {\mathrm{10}}^{-\mathrm{3}}$ mg m⁻³, respectively. Compared to the underlying layers, the CHL within the bolus is slightly elevated (∼3.6 %), while the DO is significantly higher by approximately 285 %. The thickness of the layer between these two isopycnals varies, ranging from less than 40 m, and the thickness of the trapped bolus is approximately 100 m, indicating a distinct water mass that is also associated with a low BVF. The observed elevated BVF around the bolus suggests that this is a stable water mass isolated from the surrounding water column rather than a result of vertical mixing. This lens is slightly warmer (∼22.3 °C) and more saline (∼40.43) than other waters within the same isopycnal range along the transect (Fig. 6c, d, and f). While its signature was not reflected in CHL (Fig. 6g), the bolus is also detectable in backscatter (Fig. 6h), with a concentration nearly 11 % higher than the surrounding waters (Fig. 6h). This bolus is likely outflow water from the Gulf of Aqaba, which might be advected into the region by the southward flow and subsequently captured and recirculated by the observed AE (Fig. 6a). Only a few studies are available on the water mass characteristics of the Gulf of Aqaba (Manasrah, 2002; Manasrah et al., 2004), suggesting that the upper 300 m of the gulf exhibit conditions similar to those found in the upper 100 m of the NRS during winter, with temperatures ranging from 20.4 to 22.4 °C for the salinity between 40.3 and 40.7.

The near-surface temperature continued to decrease through March, while salinity increased within the surface layer, reaching a maximum salinity of 40.4 and a density anomaly of about 28.4 kg m⁻³ in late March (Fig. 4). The corresponding glider section of temperature and salinity for this period (Fig. 7c and d), 26–29 March, shows that the coolest, saltiest, and densest water occurred in the center of a cyclonic eddy about 40–45 km offshore, where the dense isopycnals (>28.2 kg m⁻³) outcrop at the surface. Across this transect, the near-surface temperature was less than 22.8 °C, the minimum salinity was more than 40.3, the minimum near-surface density was >28 kg m⁻³, and the minimum stratification was $<\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ (Fig. 7f). Prior to this, the shallowest depth at which the 28.2 kg m⁻³ isopycnal was observed was in the early period between 5 and 9 February at a depth of about 90 m at the offshore end of the transect. In the same transect, the isopycnal descended below 250 m nearshore. The mixed layer depth along this transect ranged from as shallow as 12 m to as deep as 148 m at the inshore end of the transect, on the periphery of the eddy circulation. In the eddy center, the mixed layer extended to the depth of the 28.3 kg m⁻³ isopycnal, where stratification due to the uplifted pycnocline impeded deeper mixing. The isopycnal uplifting appears to be in the center of a cyclonic feature, where the geostrophic velocity (Fig. 7a and b) is northward nearshore (0.5 m s⁻¹) and southward offshore (0.2 m s⁻¹).

Figure 7 Same as for Fig. 5 but for the period 26–29 March 2019.

As in the earlier sections, small-scale structures are apparent and potentially important to biogeochemical processes. At about 20 km offshore, a low-DO and low-CHL feature extends upward from about 200 m to at least 100 m in depth (Fig. 7e–g). The low-DO concentration extends to shallower depths. The uplift of isopycnals affects the biogeochemical processes by bringing low-DO and low-CHL waters into the Z_eu. This process modulates nutrient, carbon, and DO availability and ultimately affects primary production. Phytoplankton growth depends on the nutrients and light availability. The low-CHL waters typically indicate nutrient-depleted conditions at the surface, while the low-DO waters in the deeper layers are generally enriched with remineralized nutrients such as nitrate, phosphate, and silicate (Garcia et al., 2019b). In this case, the low-CHL and low-DO waters have reached waters have reached ∼60 m, penetrating the Z_eu, which as reported by Zarokanellos and Jones (2021). When these nutrient-available these nutrient-available waters reach the Z_eu, they can stimulate phytoplankton blooms, enhancing primary production (Falkowski et al., 1998). The uplift of the 28.3 isopycnal (∼60 m) due to the presence of the cyclonic eddy (Fig. 7) also influences the nutrient availability (Zarokanellos and Jones, 2021; Kürten et al., 2019). This mesoscale eddy activity in the region often drives the shift in the phytoplankton community (Kürten et al., 2019).

The presence of the CE leads to an uplift of the isopycnals about 50 km offshore (Fig. 7). On the shoreward periphery of this eddy (∼20 km offshore), both CHL and DO penetrate to a depth of up to 250 m, well below the mixed layer and the Z_eu depth (Fig. 7g and e). This subduction occurs near the nearshore reversal of the flow. Unfortunately, no backscatter data were acquired during this period (Fig. 7h).

Immediately following this cool period, warmer, fresher water began to appear in the nearshore region of the glider line (Fig. 8c and d). Shallowing of lower-oxygen and hence nutrient-rich water offshore between 35 and 50 km was observed in late March and subsided with the weakening of the cyclonic circulation (Fig. 8e). The cyclonic circulation remained (Fig. 8a), but it entrained the warmer, fresher water from the south, as is evident in the composite SST image from 30 March to 6 April (Fig. 3k). This warming and freshening continued through the remainder of the glider deployment (Fig. 4a). Freshening was evident in the nearshore half of the glider line, where salinities fell below 40.3, but they remained higher along the offshore half of the line. Denser water subsided, except in the outer half of the transect. Although the direction of the geostrophic velocity was similar in pattern to the previous period (Fig. 7e), the magnitude of the nearshore flow intensified by 0.2 m s⁻¹, while the offshore flow was similar to the previous transect (Fig. 8b). Figure 8g shows that the CHL concentration between 35 and 50 km offshore decreased as warmer water was entrained near the coast.

Subduction is also present between 30 and 10 km onshore and is clearly observed in the BVF panel and in CHL and DO (Fig. 8b, e, and g). The subduction displaces the pigments of CHL and DO deeper than 200 m nearshore; in the offshore area, these two parameters are distributed well in the first 100 m.

Figure 8 Same as Fig. 5 but for the period 29 March–2 April 2019.

4 Discussion

The NRS has a dynamic and complex three-dimensional circulation with significant seasonal variability influenced by strong atmospheric forcing through wind stress and air–sea buoyancy fluxes. Direct observations and modeling experiments have both captured the formation of locally produced intermediate (RSOW) and deep water (RSDW) masses and their interactions with the adjacent gulfs of Suez and Aqaba (Table 2; Asfahani et al., 2020; Sofianos and Johns, 2003; Papadopoulos et al., 2015). Two main thermohaline cells are associated with water mass formation and influenced by mesoscale dynamics, wintertime cooling, and deep convection. Numerical simulation studies suggest that the cyclonic gyre is the most probable site of RSOW formation (Yao et al., 2014a; Sofianos and Johns, 2003).

Table 2 Summary of the major conclusions from the related studies relative to the formation of the RSOW in the NRS.

Unlike previous observations and interpretations of the NRS (e.g., Asfahani et al., 2020; Papadopoulos et al., 2015; Yao and Hoteit, 2018), this study observed a reversal of the currents in the eastern half of the basin that prevented the inflow of warmer, fresher water from the south. During this phase, the upper layer of the NRS became relatively homogeneous, and near-surface water along the glider line reached its maximum salinities and densities (Figs. 5d and 7d).

Figure 9 The first mode of the EOF based on the SLA weekly mean in 2019. (a) The spatial pattern of the first mode. (b) Time series graph of the first mode for four (4) subsequent years: 2016 (black), 2017 (red), 2018 (green), and 2019 (blue) from January to May.

To evaluate the similarities to and differences from previous years, an EOF analysis was performed of the SLA data between 26 and 28° N over a period of 4 years that was considered efficient (from 2016 to 2019; Fig. 9). The first mode of the EOF describes 86.6 % of the SLA variation (Fig. 9a). In the years 2016–2018, the EOF of the SLA showed a relatively positive or neutral pattern during the winter–spring transition period, which continued until May, when the EOFs typically decreased (Fig. 9b). This late spring decrease generally coincides with the transition from a net negative air–sea heat flux to a net positive flux (see Fig. 2b). For 2019, the first mode of the EOF showed a distinct increase in late January through mid-February and then became negative through mid-March, in contrast to the pattern in previous years. This negative phase is consistent with the period when the circulation was anticyclonic (Fig. 9b). The flow of warmer, fresher water from the south was apparently blocked during this period, and the temperature became relatively homogeneous in the NRS, as was previously observed in the CRS (Zarokanellos et al., 2017). Recently, a study by Mohamed and Skliris (2025) showed that the annual climatology of the sea level is generally higher on the eastern boundary of the RS compared with the western boundary, where many areas with isolated patches of higher or lower values of the SLA indicate mesoscale activity. The maximum values correspond to regions where AEs are present. In our study, the negative phase of the EOF analysis (Fig. 9b) aligns with the presence of AE in the NRS.

We examined the relationship between the first mode of the EOF and atmospheric forcing using a time-lagged correlation. No clear correlation was discernible from this analysis (not shown). However, in comparison with the previously published observations from 2016 (Asfahani et al., 2020), the period of negative average heat flux lasted longer into the spring. Thus, there is an overall difference in the duration of the negative heat flux between the 2016 observations and the observations in 2019.

Figure 10 Comparison between the observed transect-averaged MLD from the glider (black) and PWP-simulated MLD during the (1) cooling phase (red); the (2) cool, salty, and dense phase (blue); and the (3) warming–freshening phase (green). The black + symbol shows the MLD for each dive, interconnected by the black dashed line.

The surface layers responded to the heat loss with decreased temperatures and increased salinity and density. The cumulative effect of the cooling through the entire winter period resulted in the formation of the densest surface waters in late March, when the difference in temperature and salinity between the surface and the deep layers was at a minimum. Following this cooling phase, the net heat flux fluctuated around zero. Then, in early April, the weakening of the atmospheric forcing, the transition to positive heat fluxes, and the re-stratification due to advection from the south resulted in a near-surface temperature increase of 1 °C and a salinity decrease of 0.2, both of which contributed to the near-surface density decrease in April (Fig. 4). The warmer, less saline, and thus lighter water from the south spread into the area, and during the re-stratification it overrode the denser waters, isolating them from additional direct ocean–atmosphere interaction.

The water underlying the more buoyant surface water, which lies along the 28.2 kg m⁻³ isopycnal, extended from the surface in the middle of the transect to the approximately 200 m nearest the coast (Fig. 8). This recently exposed subsurface water spreads in the basin, and its signal can be detected in the CRS, as mentioned by Zarokanellos and Jones (2021). This water results from the northward advection of Gulf of Aden Surface Water (GASW), which is subjected to evaporation along its entire transit of the RS. Winter cooling throughout the entire period further modifies the transported surface water in the NRS. During this particular year, the presence of the AE in the NRS temporarily blocked advection from the south, contributing to the surface waters' extended exposure to the atmosphere. Near the end of the cooling period, when the surface water reached its maximum density, the cyclonic circulation was re-established, contributing to the inflow of buoyant GASW which overlaid the denser water.

Figure 11 Characteristics of the subducted feature in transect 2 on 4–9 February (Fig. 5). Panels (a) and (b show the concentrations of CHL and DO as a function of density for profiles between 20 and 40 km offshore. Panel (c) shows the integrated CHL as a function of distance offshore. The solid green line is the total integrated CHL between the surface and 500 m. The dashed green line is the integrated CHL between the surface and 150 m, and the blue line is the integrated CHL between 150 and 500 m. Panel (d) shows the percentage of the integrated CHL from 0 to 150 m (dashed green line) and from 150 to 500 m (dashed blue line).

In this study, the PWP model was applied to subsets of the observational period to further understand the relationship between the local heat flux and the advection of water from the south. The PWP model used daily surface heat flux and wind stress, as a pronounced diurnal cycle was not evident in the observed salinity and temperature data. While this simple one-dimensional model cannot capture the spatial variability of the water column structure or the atmospheric forcing field, it effectively illustrates the role of atmospheric forcing in driving the seasonal evolution of the mixed layer in the absence of these complexities. In contrast, Krokos et al. (2022) used a four-dimensional MIT-GCM model to investigate the spatial and seasonal evolution of mixed layer variability across the entire RS, highlighting the critical role of atmospheric forcing, particularly through its influence on mixed layer temperature. This broader modeling approach supports the atmospherically driven dynamics demonstrated by our PWP findings.

Figure 10 shows the evolution of the MLD in the upper 200 m based on the PWP simulation. Three separate PWP simulations were performed, initiating each simulation at the onset of one of the three phases determined from the in situ observations. The initial temperature and salinity profiles for the model run were taken from the glider section nearest the initiation point of the run. The cooling phase extended from 1 February until 8 March (35 d). During this period, the first presences of a mesoscale cyclonic eddy and later an anticyclonic eddy in the study area were observed. In the same period, the simulated MLD constantly deepened, reaching a maximum depth of 162 m at the end of the cooling phase.

Figure 12 Conceptual diagram of RSOW subduction and its biogeochemical impact on the NRS. Topographic and bathymetric representations of the study area where physical and atmospheric processes superimposed onto the map indicate the general cyclonic circulation (grey arrow), the influence of atmospheric forcing (blue arrows), and heat fluxes (grey spirals; a). The panels on the right display the glider section of CHL (b), backscatter at 650 nm (c), and DO (d). The red rows correspond to the uplifting of the isopycnal and the yellow arrows to the subduction of the newly formed RSOW. Subducted CHL, backscatter at 650 nm, and DO are traced below the Z_eu (∼120 m) associated with carbon export and ventilation of the deeper layers. The green isopleth of 180 µmol kg⁻¹ represents the nitracline (d).

During the dense phase, the observed MLD shows large fluctuations and mismatches with the simulated MLD. The one-dimensional model failed to capture the shallowing of MLD during the dense phase. The observed discrepancy is based on the pycnocline depth that shallowed substantially, such that the dense pycnocline intersected with the surface. Lastly, the PWP-simulated MLD during the warming phase also shows a discrepancy with the observed MLD. During this phase, the observed MLD rapidly shallowed to less than 45 m, and the mean MLD was about 20 m. The shallowing of the MLD during the warming and freshening phases was attributed, in part, to the northward flux of the buoyant Gulf of Aden Water. But the shallowing of the pycnocline due to the cyclonic eddy also contributed to the shallowing of the MLD, as seen clearly in Fig. 8. In addition, the shallowing of the MLD due to the presence of the CE facilitates exporting the surface water masses to the deeper layers below the euphotic zone, as indicated from the backscatter (Fig. 5h).

Water mass subduction along isopycnals is a component of the formation of RSOW during winter and a contributor to the carbon flux from the Z_eu to the interior of the RS. As shown in Fig. 5, a water mass containing elevated CHL and DO can extend well below the mixed layer and Z_eu, although only one example with this feature (5–9 February, Fig. 5) is shown in this paper. Subducted water was evident in the glider deployment from its deployment on 30 January through the fourth transect that was completed on 18 February. Key characteristics of this feature were the elevated CHL and DO in the 28.2 kg m⁻³ isopycnal that extended to as deep as 250 m. In Fig. 5, this feature was evident from about 45 km offshore and shoreward. Figure 11 shows the relationship between these variables and densities between 20 and 40 km offshore. A clear peak in both variables aligns with the 28.2 kg m⁻³ isopycnal. Inshore, at 45 km, the region between 150 and 500 m contains a measurable fraction of the total CHL. In addition, Fig. 11d shows that the subduction event contributes 20 % to the integrated CHL between 150 and 500 m.

Given the limitations of our observations, constrained by the Exclusive Economic Zone (EEZ) boundary, the full mechanism of the formation of this subducted layer is unclear. One possible mechanism is that either vertical mixing or sinking of particles in the western half of the NRS (Kheireddine et al., 2020) could create this feature, which is then entrained into the cyclonic circulation in this region and transported from the western side of the basin to the eastern side. Figure 12 shows a conceptual diagram of RSOW subduction and its biogeochemical impact on the NRS based on the existing observations. Regardless of the details of the mechanism, subduction is a process that needs to be considered in the physical and biogeochemical dynamics of the NRS.

5 Conclusion

The primary objective of this study was to understand the mechanisms contributing to the water mass formation of RSOW in the NRS and the associated biogeochemical responses. Our findings demonstrate that the RSOW formation is closely linked to the presence of a cyclonic eddy and intense winter cooling, consistent with previous studies (Asfahani et al., 2020; Sofianos and Johns, 2003; Papadopoulos et al., 2015; Yao et al., 2014b). A novel finding of this study is the role of water mass subduction, which, although not previously discussed, contributes not only to the RSOW formation but also to the carbon export as CHL and backscatter indicated that water has been subducted below Z_eu (Fig. 12). The observations also indicate that subduction events can significantly contribute to the ventilation of intermediate and deeper waters, thereby affecting the overall oxygen budget of the RS. During these events, oxygen-rich surface waters are transported into subsurface layers along the 28.2 kg m⁻³ isopycnal, facilitating oxygen redistribution at depth. The study also identified a transition from negative to positive heat flux and the re-establishment of northward flow along the eastern RS coast, signalling the cessation of the RSOW formation as less dense water from the south caps the denser northern waters. The presence of the AE south of the study area prevented the advection of more buoyant surface water, although we could not determine the mechanism for this reversal of CE to AE from our observations and one-dimensional model simulations. This likely influenced the observed water mass originating from the northern gulfs within the RSOW density range, consistent with Papadopoulos et al. (2015). This study highlights that multiple factors can contribute to RSOW formation (Table 2), including dense isopycnal surfacing from cyclonic circulation, vertical mixing, dense outflow from the northern gulfs, water mass subduction, the extension of dense isopycnal exposure due to blockage of buoyant flux from the south, and the eventual termination of these processes as the buoyancy flux is restored. In addition, it is clear that the submesoscale features are present in the region and contribute to the overall physical and biogeochemical dynamics of the region. To comprehensively capture the spatial and temporal dynamics of RSOW formation, future research should prioritize detailed observational and modeling studies, integrating autonomous platforms and ship-based sampling across the entire NRS basin. Such an approach would resolve the three-dimensional variability and provide valuable insights into the sources and sinks involved in RSOW formation and its biogeochemical impacts.