Mechanisms driving mesoscale latent heat flux variations and mixed layer heat content evaluation in the Northwest Tropical Atlantic

In this study, a high-resolution ocean-atmosphere coupled simulation is used to assess the effects of sea surface temperature (SST), surface currents, and ocean vertical stratification on the spatial variability of latent heat flux (LHF) and the stability of the marine atmospheric boundary layer (MABL) in the Northwest Tropical Atlantic during January and February 2020. The analysis focuses on the ocean mesoscale (O(50–250 km)) across the Northwest Tropical Atlantic (referred to as the EURECA region in this study) and within three sub-regions characterized by different ocean dynamical regimes: Amazon, Downstream, and Tradewind. Results indicate that the coupling between SST and wind speed (and specific humidity) is stronger (weaker) in the Amazon and Downstream regions, influenced by the warm coastal North Brazil Current eddy corridor and the Amazon River plume, than in the Tradewind region, representative of the open ocean, consistent with previous remote sensing studies. Overall, warmer SSTs are associated with increased wind speeds and variations in specific humidity, deviating from Clausius–Clapeyron expectations. We interpret this as the result of active ocean processes modifying the near-surface atmosphere, enhancing vertical motion in the MABL, and transporting momentum and drier air from the free troposphere toward the surface. To further investigate the impact of mesoscale SST features on LHF, we apply a linear, SST-based downscaling method. Results show that these mesoscale SST structures induce a substantial increase in LHF, 46.8 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ on average in the Amazon and Downstream regions (warm eddy corridor). In the Tradewind region, the LHF sensitivity to SST is smaller, at about 35 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$. For the Amazon region, of the 46.7 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ change in LHF associated with SST, approximately 7.8 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ is attributed to direct mesoscale SST changes (thermodynamic contribution), while the remainder is linked to mesoscale SST-induced modifications in near-surface atmospheric circulation (dynamic contribution), mainly due to the mesoscale SST-induced humidity undersaturation imbalances. The influence of surface currents on LHF is weaker, with deviations not exceeding 15 W m⁻². Finally, we focus on the SST mesoscale anomalies linked to the Amazon freshwater plume. We find them to be persistent throughout the period of study affecting LHF by the mechanisms described above. Lateral advection and heat loss to the atmosphere tend to dilute them with their environment by the end of the period of study. This work underscores the importance of a regionalized approach to mesoscale air-sea interaction studies in the Northwest Tropical Atlantic, as LHF sensitivity to SST and surface currents exhibits strong spatial variability driven by distinct oceanic dynamics. Submesoscale LHF sensitivity to SST and currents is not addressed here and will be the subject of future research.

1 Introduction

Turbulent heat fluxes (THFs, comprised of latent and sensible heat fluxes) are related to temperature (sensible) and moisture undersaturation (latent) imbalances at the air-sea interface. When examining air-sea interactions through THFs, it is common in the literature to distinguish between the ocean's large-scale and fine-scale processes, the latter including the mesoscale (O(50–250 km)) and submesoscale (<O(50 km)) components. These fine-scale interactions with the atmosphere have been shown to differ significantly from large-scale processes (Chelton and Xie, 2010; Small et al., 2019; Gentemann et al., 2020; Conejero et al., 2024). At large scales, atmospheric dynamics predominantly drive ocean variability (Gill, 1982). However, at scales smaller than about 250 km, the ocean actively influences the near-surface atmosphere, affecting air temperature, frictional stress, and the stability of the marine atmospheric boundary layer (MABLH; Small et al., 2008). Among THFs, this study focuses on latent heat flux (LHF) in the Northwest Tropical Atlantic, as it provides a direct link between atmospheric dynamics and thermodynamics. Indeed, the process of seawater evaporation cools the ocean surface, while the heat released during the subsequent moisture condensation warms the atmosphere. Warm and moist air in the atmosphere can become buoyant, triggering atmospheric convection and the formation of storms.

The effects of fine-scale sea-surface temperature (SST) variability on the near-surface atmosphere and air-sea heat fluxes, have been investigated using satellite products (Bishop et al., 2017; Fernández et al., 2023), in-situ measurements (Acquistapace et al., 2022; Iyer et al., 2022; Fernández et al., 2024), atmospheric models (Borgnino et al., 2025) and coupled models (Small et al., 2019). In the literature, two primary mechanisms of lower atmospheric response to SST features have been identified: the downward momentum mixing (DMM; Hayes et al., 1989; Wallace et al., 1989) and the pressure adjustment (PA; Lindzen and Nigam, 1987). The DMM mechanism consists in the destabilization of the MABL over a warm SST anomaly, which enhances vertical mixing (Fig. 1a). This process facilitates the entrainment of drier air from the free troposphere into the MABL, thereby increasing surface winds and intensifying LHF (Acquistapace et al., 2022; Borgnino et al., 2025). Conversely, cold SST anomalies suppress vertical mixing, leading to reduced surface winds and lower LHF. Thus, DMM provides a top-down mechanism by which fine-scale SST variability influences the near-surface atmosphere, a process known as thermal feedback (TFB; Renault et al., 2019b, 2023). PA, on the other hand, predicts that surface wind convergence (divergence) occurs over SST maxima (minima) as warm (cold) SST cores generate local sea level pressure lows (highs). This leads to weaker winds over SST extrema, resulting in lower LHF (Pasquero et al., 2021). The influence of these mechanisms has been observed across various regions of the World Ocean on timescales ranging from hours to days and months (Small et al., 2023). Foussard et al. (2019) found that PA tends to dominate where surface winds are well-coupled to upper-level winds, while DMM prevails in neutrally stable lower tropospheric conditions where SST effectively modulates surface winds (Desbiolles et al., 2023). Warm mesoscale eddies and fronts have been shown to enhance LHF via DMM in several regions, including the Gulf Stream (Minobe et al., 2008), the Kuroshio Extension (Xu et al., 2011; Ma et al., 2015; Chen et al., 2017), the South China Sea (Liu et al., 2018, 2020), and the Agulhas (O'Neill et al., 2005) and Malvinas (Villas Bôas et al., 2015; Leyba et al., 2017) currents. Meanwhile, PA has been shown to impact cloud and precipitation patterns in the cold wake of tropical cyclones (Ma et al., 2020) through a cross-track secondary circulation (Pasquero et al., 2021). Finally, Conejero et al. (2025) showed that both PA and DMM are diminished when the model grid spacing is reduced and the full range of mesoscale structures is resolved due to a stronger submesoscale-induced atmospheric frontogenesis. In the Northwest Tropical Atlantic, Fernández et al. (2023) found that DMM dominates over PA using multiple satellite products.

Figure 1 (a) Schematic representation of the Downward Momentum Mixing (DMM) mechanism, adapted from Meroni et al. (2020). (b) Schematic representation of the Current Feedback (CFB) mechanism. Details on how to compute the CFB-induced wind anomaly (ΔU_CFB) are provided in Sect. 3.1

In addition to SST effects, surface currents, also influence wind stress and thus surface winds and LHF via a bottom-up process known as current feedback (CFB; Bye, 1985; Chelton et al., 2001; Renault et al., 2016a, b, 2019a). Here, we focus on the CFB-induced surface wind response and its impact on LHF. When surface currents (black arrow in Fig. 1b) and winds (blue arrow) align (left side of the eddy in Fig. 1b), surface stress becomes reduced and an anomalous surface wind develops (orange arrow), reinforcing the prevailing wind field (Renault et al., 2016b). This enhances relative wind speed (the difference between surface wind and surface current velocities, represented by a green arrow) which becomes larger than if we had not accounted for CFB (purple arrow) and increases LHF. Conversely, when surface currents oppose surface winds, the CFB-induced wind anomaly weakens surface winds, which could lead to smaller LHF (right side of the eddy in Fig. 1b).

The significance of CFB in eddy dynamics has been highlighted in studies using remote sensing data (Gaube et al., 2015) and high-resolution coupled simulations (Eden and Dietze, 2009; Renault et al., 2016b). Ignoring atmospheric adjustments to CFB might lead to an overestimation of eddy attenuation timescales and an underestimation of eddy amplitude and azimuthal speed. CFB might also shorten eddy lifetimes. Indeed, the observed composite life cycle reconstructed from satellite altimetry in Renault et al. (2016b), consists of a rapid early intensification, a prolonged slow CFB-induced decay, and an abrupt collapse at the end. This would imply that CFB systematically acts as an eddy-killer, transferring energy from mesoscale eddies to the atmosphere (Bye, 1985; Dewar and Flierl, 1987; Anderson et al., 2011; Oerder et al., 2018; Renault et al., 2016b). Furthermore, an inaccurate representation of CFB in numerical models has implications for biogeochemistry. In oligotrophic regions, mesoscale processes enhance the upward transport of limiting nutrients, supporting biological production (Martin and Richards, 2001; Gaube et al., 2013). In coastal upwelling systems, such as along the Californian coast, eddies modulate biological productivity by subducting nutrients below the euphotic zone and advecting biogeochemical material offshore (Gruber et al., 2011; Nagai et al., 2015; Renault et al., 2016a). Misrepresenting CFB may lead to an overestimation of nutrient quenching and offshore transport, thereby impacting marine ecosystems by altering eddy amplitude, lifetime, and spatial extent.

As stated above, this paper is focused on the Northwest Tropical Atlantic. The ocean circulation in this region, particularly near the Amazon River estuary, is dominated by the North Brazil Current (NBC), as illustrated in Fig. 2a. The NBC is a strong western boundary current originating in the Equatorial and South Atlantic. Around 8° N, it separates from the coast, forming the NBC retroflection, which feeds the North Equatorial Countercurrent. The NBC system is closely linked to two major processes: the Amazon freshwater advection in the open ocean (Reverdin et al., 2021) and the periodic formation of mesoscale eddies known as NBC rings (Johns et al., 1990; Richardson et al., 1994). These processes are interconnected, as NBC rings facilitate the offshore transport and the lateral spreading of the Amazon river plume induced by mesoscale advection and smaller-scale instabilities (Reverdin et al., 2021; Olivier et al., 2022; Coadou-Chaventon et al., 2024) before coalescing and dissipating due to their interaction with the complex topography near the Lesser Antilles (Fig. 2b; Subirade et al., 2023). This results in strong spatial heterogeneity in sea-surface salinity (SSS, Fig. 2a), which also influences upper-ocean temperature by modulating stratification. When low salinity dominates the upper-ocean stratification, the ocean mixed layer (ML), the surface layer in direct contact with the atmosphere, becomes shallower. This can lead to the formation of barrier layers (BLs), which promote temperature inversions (Mignot et al., 2012; Mahadevan et al., 2016; Krishnamohan et al., 2019; Coadou-Chaventon et al., 2024). In the presence of a BL, heat and momentum fluxes remain confined to the shallow ML, which responds more rapidly to atmospheric forcing. Consequently, the ML cools more quickly in winter and warms more rapidly in summer due to the inhibited exchange with deeper ocean layers (Miller, 1976; Sprintall and Tomczak, 1992). This leads to negative SST anomalies relative to surrounding waters over the Amazon plume and, in turn, reduced LHF in boreal winter. The opposite holds in boreal summer. However, the magnitude of this response remains debated. While observational studies suggest that BLs have a strong impact on SST (Pailler et al., 1999; Foltz and McPhaden, 2009), numerical models often fail to reproduce this effect (Breugem et al., 2008; Balaguru et al., 2012; Hernandez et al., 2016).

Figure 2 (a) Major dynamical features of the western equatorial Atlantic (arrows) overlaid on the averaged sea-surface salinity (SSS) field from 17–19 February 2020. The SSS field is derived from the SMAP-SSS Level 3, version 4.0, 8 d running-mean gridded product (Boutin et al., 2021). The black box delineates the EURECA region, which corresponds to the simulation domain used in this study. (b) Seafloor depth from the simulation used in this study which is derived from the GEBCO 2020 global bathymetric dataset (Conejero et al., 2024). The three boxes delineate the three sub-regions within the EURECA domain assessed in this paper: Amazon, Downstream and Tradewind and the black contour represents the 100 m isobath. Numbers 1 and 2 mark the locations of Trinidad and Tobago and a region near Barbados, respectively. These geographical references, along with the Lesser Antilles marked in panel (a), are used throughout the main text. Detailed descriptions of the simulation used in this study and of the sub-regions follow in Sects. 2 and 3 respectively.

The objective of this study is to investigate how these processes (TFB and CFB) together with the Amazon plume affect LHF variations in the Northwest Tropical Atlantic, complementing previous research based on satellite/reanalysis products (Fernández et al., 2023) and in-situ observations (Fernández et al., 2024). To achieve this, we utilize the high-resolution regional coupled simulation described in Conejero et al. (2024), fully resolving the ocean mesoscale. Such a simulation provides a view in three dimensions of both the ocean and the atmosphere, allowing for statistical robust results. This contrasts with satellite observations, which only provide ocean or atmosphere surface fields. In-situ observations from vertical transects, in turn, do provide high-resolution information about vertical structures, but they are sparse in time and space.

This research is conducted within the framework of the Elucidating the Role of Cloud-Circulation Coupling in Climate – Ocean Atmosphere (EUREC⁴A-OA, http://www.eurec4a.eu, last access: 12 January 2026) and the Atlantic Tradewind Ocean-Atmosphere Mesoscale Interaction Campaign (ATOMIC, https://psl.noaa.gov/atomic/, last access: 12 January 2026) field experiments. The paper is structured as follows: Sect. 2 describes the simulation configuration. Section 3 presents the methods used to analyze model data. The main results are discussed in Sect. 4, followed by conclusions and future perspectives in Sect. 5.

2 Data

This study utilizes the EURECA ocean-atmosphere coupled simulation (Conejero et al., 2024). The ocean component is based on the Coastal and Regional Ocean COmmunity (CROCO) model (Shchepetkin and McWilliams, 2005; Debreu et al., 2012), while the atmospheric component employs the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). The two models are coupled via OASIS (Craig et al., 2017), which performs the grid interpolation and temporal averaging for property exchanges between the two model components every hour. Specifically, the ocean model provides SST and surface currents to the atmosphere, while WRF returns surface heat, momentum, and water fluxes to CROCO. Further details on the EURECA simulation configuration are available in Conejero et al. (2024).

The EURECA simulation spans from January 2019 to June 2020, though this study focuses on the January–February (JF) 2020 period. The CROCO domain extends from 5.5 to 15.5° N and 62 to 52° W, with a horizontal resolution of 1 km. Initial and lateral open boundary conditions are provided by the “Antilles” simulation, which employs the same coupled configuration but at a coarser resolution over a larger domain, including parts of the Caribbean Sea (Conejero et al., 2025).

The WRF domain in EURECA is slightly larger than the ocean domain to mitigate sponge effects. The atmospheric model runs at a horizontal resolution of approximately 2 km, with outputs available on 40 η vertical levels. To analyze the first 2000 m of the atmosphere, these levels are linearly interpolated to uniform 100 m vertical spacing. Initial and boundary conditions, provided by the “Antilles” simulation, are updated every 3 h. Bulk formulations (Fairall et al., 2003) are used to estimate freshwater and turbulent fluxes, which are then fed into the ocean model. The CFB effect is implemented in both the surface and planetary boundary layer schemes, following Renault et al. (2019a). WRF variables at multiple vertical levels are stored every 3 h, while surface variables are recorded hourly.

3 Methodology

Following Fernández et al. (2024), we remove the diurnal cycle from all variables before analyzing LHF sensitivity to the surface ocean. This is achieved by computing daily means from the WRF/CROCO outputs, and all subsequent calculations are performed on these averaged variables. It is important to note that this procedure filters out part of the ocean submesoscale variability. Consequently, this paper deals exclusively with the ocean mesoscale. We briefly return to this point in Sect. 5.

All analyses are performed in four regions. First, we consider the full simulation domain, referred to as the EURECA domain (5.5–15° N, 62–52° W). Additionally, we examine three distinct sub-regions: Amazon (7–9° N, 56–53° W), Downstream (9–11° N, 59–56° W) and Tradewind (12–14° N, 58–52.5° W). Figure 2b displays three boxes delineating the sub-regions. The three subdomains are characterized by different ocean dynamics and air-sea interactions. The Tradewind sub-region, representative of the open ocean, is relatively quiescent, whereas the Downstream sub-region, closer to the coast, exhibits enhanced fine-scale ocean activity. In turn, the warm and fresh surface waters associated with the Amazon plume are advected into the Amazon sub-region as shown in following sections. These regional differences have been highlighted in previous studies based on remote sensing and reanalysis data (Fernández et al., 2023), as well as in-situ observations (Fernández et al., 2024). To ensure consistency with these studies we focus on the JF 2020 period. In some cases, the results are only presented for February 2020 to highlight the impacts of the Amazon freshwater plume, which only reaches Amazon by mid-February. This is clearly indicated.

3.1 Coupling Coefficients

Following Renault et al. (2016b), Renault et al. (2019b), and Conejero et al. (2024), we estimate several air-sea coupling coefficients as the statistically significant slope (determined via a two-sided t test) of the linear regression between the binned distributions of mesoscale anomalies from two variables. For all coefficient calculations, we exclude mesoscale anomalies located on the continental shelf (regions where the seafloor is shallower than 100 m, see Fig. 2b) and the Lesser Antilles (west of 60.25° W), as these areas can be affected by orography, land-sea interactions, and coastline effects (Desbiolles et al., 2014).

To compute mesoscale anomalies we use a combination of time and spatial filters. In order to remove weather-induced synoptic variability from the atmospheric variables (winds, specific humidity, and air temperature), we first apply a 29 d running mean as in Chelton et al. (2007), Renault et al. (2019b), and Conejero et al. (2024). To isolate the mesoscale band, we apply a band-pass isotropic Gaussian filter with a 50–250 km cutoff length (Renault et al., 2019b), to the 29 d running mean dataset and keep the scales between 50–250 km. After performing the filtering for the whole EURECA domain, we select the anomalies in Amazon, Downstream and Tradewind to operate on them separately.

To characterize the effect of surface currents on near-surface winds (CFB) we compute s_w, the coupling coefficient between surface current vorticity and surface wind curl anomalies (Renault et al., 2016b, 2019b). Note that, since mesoscale currents are nearly in geostrophic balance (and therefore non-divergent), s_w effectively isolates the CFB from the TFB at the mesoscale (Renault et al., 2019b). Using s_w, the CFB-induced wind anomaly (ΔU_CFB) reads:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\mathit{U}}_{\text{CFB}}=s_w{\mathit{U}}_{\mathrm{o}},\end{array}$$

where U_o stands for surface currents. Recall that ΔU_CFB is represented in orange arrows in Fig. 1b.

Finally, we also calculate the coupling coefficients of near-surface wind speed (s_u), near-surface specific humidity (s_q) and near-surface atmospheric temperature (s_t) mesoscale anomalies with respect to SST mesoscale anomalies. Table 1 provides a summary of the coupling coefficients described above.

Table 1 Overview of the coupling coefficients.

3.2 LHF Sensitivity to SST and Surface Currents

To evaluate how LHF responds to SST and surface currents, we compute multiple LHF datasets using the COARE3.5 algorithm (Edson et al., 2013). First, we calculate LHF_U, which represents the LHF dataset obtained using the atmospheric variables of the first WRF vertical level (surface winds, air temperature and specific humidity) together with SST. Note that we do not consider relative winds (i.e. the difference between surface winds and surface currents) to compute LHF_U. Additionally, we compute LHF_LR, which corresponds to LHF computed with the same variables as LHF_U, but smoothed with a Gaussian low-pass filter (cutoff length of 250 km). Again, we do not consider relative winds in the calculation of LHF_LR either.

We then apply the LHF downscaling algorithm developed by Fernández et al. (2023). Given a smoothed variable (Ψ_LR), we reconstruct a new dataset incorporating the finer-scale SST features (Ψ_HR) as:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathrm{\Psi }}_{\text{HR}}={\mathrm{\Psi }}_{\text{LR}}+s_{\mathit{\psi }}\mathrm{\Delta }\text{SST}.\end{array}$$

Here, ΔSST represents the SST correction, which accounts for deviations of the high-resolution SST field from the coarse smoothed SSTs. To ensure that the domain-averaged SST correction is zero and that the area-weighted means of the variables remain conserved, ΔSST is computed as:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }\text{SST}=\left(\text{SST}-\stackrel{\mathrm{‾}}{\text{SST}}\right)-\left({\text{SST}}_{\text{LR}}-\stackrel{\mathrm{‾}}{{\text{SST}}_{\text{LR}}}\right).\end{array}$$

$\stackrel{\mathrm{‾}}{\text{Overbars}}$ denote the spatial average over the region of study. Thus, to compute LHF_HR, we statistically downscale each of the variables driving LHF (surface wind, specific humidity air temperature and SST) using its corresponding coupling coefficient, as detailed above, and we then apply COARE3.5 with the downscaled variables. Note that we do not use relative winds for LHF_HR either.

One shortcoming of this downscaling algorithm is the assumption that the reconstruction of the finer-scale fields exclusively depends on SST changes, which is not always the case (i.e. specific humidity depends on wind speed; Small et al., 2019). However, the algorithm has been shown to improve LHF estimations by a factor of two in an ensemble of SST-forced WRF atmospheric simulations (Fernández et al., 2023) and allows to isolate the contributions from the different LHF controlling variables (i.e. surface winds, specific humidity etc) when the ocean SST mesoscale effect in them is considered. These facts along with its simplicity of implementation encourage us to use it here despite this shortcoming.

To isolate the thermodynamic contribution (i.e. LHF variations solely due to SST changes via modifications of the saturation specific humidity maintaining relative humidity at its smoothed value) to LHF sensitivity, we compute an additional LHF dataset, denoted as LHF_therm. In this dataset, air temperature and SST are modified adding ΔSST (no s_t involved here in order to keep the air-sea temperature imbalance constant, although a proper downscaling of air temperature using Eq. (2) does not significantly modify the results, suggesting a weak role of air temperature in LHF sensitivity to SST) and wind speed and relative humidity remain with their smoothed values. The specific humidity required to obtain LHF_therm is derived from the smoothed relative humidity and the specific humidity of saturation computed with SST_HR via the Clausius–Clapeyron equation. This ensures that the specific humidity variations for this LHF subset are only due to SST mesoscale changes themselves without the effects of the atmospheric-induced modifications of mesoscale SST structures (i.e. entrainment of drier air from the free troposphere associated with DMM).

Moreover, we compute LHF_therm−U to further distinguish, within the dynamic contribution, the effects of wind speed and relative humidity. The dynamic contribution represents LHF changes associated with the mesoscale SST-induced modifications of the near-surface atmosphere (winds and relative humidity, the contribution of air-temperature is negligible and not assessed here). Thus, in LHF_therm−U, air temperature, and surface winds are downscaled, while specific humidity is obtained as in LHF_therm (i.e. using the smoothed relative humidity and the downscaled saturation specific humidity). This means that the only change between LHF_HR and LHF_therm−U is that specific humidity in the latter is scaled maintaining relative humidity at its smoothed value whereas in the former, specific humidity is downscaled with Eq. (2). Thus, LHF_therm−U accounts for both the thermodynamic contribution to LHF sensitivity to SST and the effect of SST-induced mesoscale surface wind changes, but it does not account for any effect of SST anomalies on relative humidity.

Finally, to assess the contributions of surface currents in LHF variations, we consider two additional LHF datasets. LHF_orig is computed using the original WRF variables in COARE3.5, but replacing the first vertical model-level wind speed with relative winds. Thus, the difference LHF_orig−LHF_U represents the variations in LHF associated with the consideration of relative winds instead of just surface winds when computing LHF. Meanwhile, LHF_no-CFB stands for LHF computed with relative winds, but with the CFB-induced wind anomaly (ΔU_CFB, Eq. 1, represented in orange arrows in Fig. 1b) removed component-wise. Hence, LHF_orig−LHF_no-CFB contains the CFB effect on LHF. Table 2 contains a summary of the different processes targeted and the LHF dataset differences used to isolate them.

Table 2 Processes targeted and corresponding LHF datasets used to isolate them.

3.3 Mixed Layer Heat Budget

As shown later, the Amazon sub-region is crossed by the fresher Amazon River plume waters, which affect ocean stratification and the SST mesoscale anomaly field in the sub-region. To provide a broader insight on the linkages between SST anomalies and the processes leading to heat redistribution in the Amazon plume, we present a mixed layer heat budget analysis in the Amazon sub-region. Based on observational studies in the region (Reverdin et al., 2021), we set the boundary of the Amazon plume at the 35 psu isoline (Amazon plume waters <35 psu).

The temperature equation, vertically integrated down to the mixed layer depth reads (Vialard and Delecluse, 1998):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\text{Total tendency}}{\underbrace{{〈{\partial }_{t}T〉}_H}}=\underset{\text{Horizontal advection}}{\underbrace{{〈-u{\partial }_x\mathrm{T}-v{\partial }_{y}T〉}_H}}-\underset{\text{Vertical advection}}{\underbrace{{〈w{\partial }_{z}T〉}_H}}\\ \text{(4)}& {\displaystyle }& {\displaystyle }+\stackrel{\text{Forcing }\mathcal{F}}{\overbrace{\underset{\text{Atmospheric forcing}}{\underbrace{{\displaystyle \frac{Q_{\mathrm{s}}\left(\mathrm{1}-F_H\right)+Q_{\text{ns}}}{{\mathit{\rho }}_{\mathrm{0}}{C}_{p}H}}}}}}+\stackrel{\text{Mixing}\mathcal{D}}{\overbrace{\underset{\text{Entrainment}}{\underbrace{{\displaystyle \frac{{\partial }_{t}H}{H}}\left(T_H-{〈T〉}_H\right)}}+\text{Residual}}}\end{array}$$

where angle brackets (〈〉) indicate integration down to the base of the mixed layer. To facilitate interpretation, the mixed layer depth (MLD) is denoted as H in Eq. (4). In this equation, T represents the ocean temperature, u, v, and w denote the zonal, meridional, and vertical currents, respectively, and T_H is the temperature at the base of the mixed layer. The MLD is computed using a density threshold criterion of $\mathrm{\Delta }\mathit{\sigma }=\mathrm{0.01}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, following Gévaudan et al. (2021). This criterion yields MLD values consistent with the in-situ observations' in Fernández et al. (2024). In Eq. (4), Q_s represents the solar component of the total heat flux, primarily shortwave radiation (SW). F_H denotes the fraction of solar radiation reaching the ML base.

Note that the flux sign convention differs when computing the MLD heat budget compared to the atmospheric convention: fluxes directed from the ocean to the atmosphere are considered negative, as they contribute to ML cooling. Finally, the EURECA simulation employs the COARE3.0 bulk formulae (Fairall et al., 2003) to compute turbulent heat fluxes (THFs). For consistency, model-derived THFs are used when calculating the ML heat budget. However, we used COARE3.5 (Edson et al., 2013) to assess LHF sensitivity to SST and surface currents. This does not lead to any inconsistency in the results, since we always compare between LHFs computed with the same algorithm. In addition, we verified that the differences between both algorithms are small compared to the LHF difference values obtained in this article.

Finally, the residual term in Eq. (4) accounts for horizontal and vertical diffusion, as well as numerical errors associated with the computation of the other terms. The reader is referred to Appendix A for further details on the most appropriate way to compute the residual to minimize numerical errors. Additionally, we explore various criteria to determine whether vertical or horizontal diffusion dominates the residual in Appendix B.

Another important quantity used to analyze the effects of water temperature in vertical stratification is the base of the isothermal layer (THERM). As in Gévaudan et al. (2021), we estimate it as the depth at which the water temperature is 0.2 °C lower than the 10 m depth level temperature. Therefore, the barrier layer thickness (BLT) results from the difference between THERM and MLD. Finally, to quantify the relative importance of salinity in ocean stratification, we use the Ocean Stratification Strength (OSS) indicator (Maes and O'Kane, 2014; Gévaudan et al., 2021):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\text{OSS}={\displaystyle \frac{N_S^{\mathrm{2}}}{N^{\mathrm{2}}}}.\end{array}$$

Here, N² represents the Brunt–Väisälä frequency, while $N_S^{\mathrm{2}}$ denotes the Brunt–Väisälä frequency computed using a constant representative temperature, allowing only for salinity variations. OSS values greater than 50 % indicate that salinity dominates over temperature in controlling ocean vertical stratification.

4 Results

4.1 Air-Sea Interface

Figure 3 presents the JF 2020 mean state of different air-sea interface variables in the Northwest Tropical Atlantic. The shading in Fig. 3a represents SST, highlighting a warm water stripe (>27.3 °C) along the South American coast, extending from the southernmost coastal point in the domain to Trinidad and Tobago. To its east, a parallel cold SST band (<25.5 °C) borders Trinidad and Tobago and extends northward towards the Lesser Antilles. Longitude-depth transects of seawater temperature, salinity and vertical velocity at different latitudes show that the sharp topography of the region where the cold SST band lies (Fig. 2b) induces the surface upwelling of cooler deeper waters (not shown). Further offshore, SST values remain relatively homogeneous (∼27 °C), although localized anomalies are observed, such as the warmer region to the east of the Amazon subdomain. Contours in Fig. 3a show first vertical level specific humidity, which is highest near the South American coast and decreases towards the northeast.

Figure 3 JF 2020 mean of (a) SST (in °C, shaded) and near-surface specific humidity (in g kg⁻¹, contours), (b) LHF (in W m⁻², shaded) and near-surface wind speed (in m s⁻¹, contours), (c) MABLH (in m, shaded) and surface winds (arrows), and (d) SSS (in psu, shaded) and surface currents (in m s⁻¹, arrows). In all panels, from north to south, the three boxes delineate the Tradewind, Downstream, and Amazon sub-regions and numbers 1 and 2 are placed over Trinidad and Tobago and a region near Barbados respectively.

Figure 3b displays near-surface wind speed (contours), which varies between 7.5–8.2 m s⁻¹. The LHF spatial pattern (shading) closely follows the SST distribution in Fig. 3a, with the highest LHF values (>180 W m⁻²) occurring over the warm coastal stripe and in the open ocean, and the lowest values (<100 W m⁻²) located over the cold SST band extending across the continental shelf. Other patches of relatively low LHF are observed offshore over cooler SSTs such as the ones west of Barbados.

To further characterize the atmospheric JF 2020 spatial pattern, Fig. 3c presents the MABLH (shading) and surface winds (arrows). Consistent with the DMM mechanism, the shallowest MABL depths (<500 m) align with the cold SST stripe, whereas deeper MABL values occur over the open ocean, where SSTs and wind speeds are higher. The dominant wind direction shifts from easterly in the open ocean to northeasterly (trades) closer to the coast. Figure 3d displays surface currents (arrows) and SSS (shading). In the southeastern domain, strong northwesterly currents (NBC) advect a low-salinity patch (SSS<35 psu), whose shape and extent are modulated by local eddy-driven circulation. Finally, near Trinidad and Tobago, a second region of enhanced surface currents exhibits a clockwise rotation, at 10° N, 58° W. This eddy, which remains nearly stationary during JF 2020, lacks a strong salinity signature at the surface. This signature is observed when studying its vertical structure (not shown).

4.2 Coupling Coefficients

To assess the relationship between SST/surface current anomalies and the near surface atmosphere at the mesoscale, we compute the coupling coefficients as defined in Sect. 3.1. Note that these coupling coefficients will be needed to statistically downscale each of the LHF controlling variables (Eq. 2) in the following sections. Since the ocean dynamics is different between sub-regions, we compute the coupling coefficients separately in EURECA, Amazon, Downstream and Tradewind. Figure 4, illustrates the value of the coupling coefficients and the binned linear regressions are displayed Appendix C.

Figure 4 Mesoscale coupling coefficients s_w (red), s_u (cyan), s_q (orange) and s_t (green). In each dot quartet, from left to right, the markers represent the corresponding coupling coefficient in the EURECA, Amazon, Downstream and Tradewind domains respectively. The error bars depict the standard error of the slope. All the coefficients displayed here are statistically significant at the 95 % level after a two-sided t test.

The first four red markers in Fig. 4 present the spatial distribution of the slope of the linear regression between the mesoscale surface current vorticity and surface wind curl anomalies. Recall that this linear regression results in the coupling coefficient named s_w, used to assess the CFB. The intensity of s_w can be interpreted as the efficiency of the partial re-energization of the ocean through the wind response to CFB (Renault et al., 2016b, 2019b). Mesoscale s_w values range from 0.22 to 0.28, with weaker coupling in the Amazon and Tradewind subdomains (0.22 and 0.24 respectively). The strongest mesoscale coupling is found in EURECA and in the Downstream subdomain, where s_w reaches 0.28 and 0.26 respectively.

In other words, these s_w values indicate that a mesoscale eddy with a velocity of 1 m s⁻¹ induces, on average, a wind speed anomaly of 0.28 m s⁻¹ in the EURECA domain, which in turn influences LHF estimations through changes in surface winds as shown in Fig. 1b. They also show the spatial variability in the strength of the coupling. The interactions are stronger within regions with nearly stationary mesoscale eddies like Downstream than in Tradewind (open ocean).

s_u is displayed in the second four cyan markers of Fig. 4. It is positive in the four domains: warm (cold) SST anomalies increase (decrease) surface winds, which is consistent with the DMM mechanism. s_u is weaker in Tradewind (open ocean, 0.11 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) than in Amazon and Downstream (warm eddy corridor, 0.29 and 0.31 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ respectively). This feature has been previously found in Fernández et al. (2023) using satellite observations. The entire EURECA region produces a s_u of 0.26 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$.

Another key variable driving LHF is near-surface specific humidity (q). As with wind speed, we use the specific humidity from the first vertical level in the WRF simulation at about 10 m. The four orange markers in Fig. 4 show the associated coupling coefficient s_q. It exhibits weak negative values over the Amazon (−0.05 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), Downstream (−0.09 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) and EURECA (−0.03 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$; in agreement with Fernández et al., 2023), while a stronger positive value appears associated with the Tradewind subdomain (0.24 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). The underlying physical explanation is as follows: at scales larger than the mesoscale, evaporation is sufficient for specific humidity to adjust to SST according to the Clausius–Clapeyron equation. This adjustment rate was estimated to be 1.3 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ through linearization of the Clausius–Clapeyron equation in Fernández et al. (2023). However, at the mesoscale this equilibrium does not hold, as q does not have sufficient time to adjust to SST variations, resulting in a weak correlation between the two variables. For example, in Amazon, the mesoscale s_q is negative meaning that specific humidity decreases with increasing SST. This behavior may reflect the entrainment of colder and drier air from the free troposphere due to the DMM mechanism. On the contrary, even though s_q is positive in Tradewind, when examining in detail the binned distribution, we find a negative (positive) slope for negative (positive) SST mesoscale anomalies separately (Fig. C3d of Appendix C). This fact highlighs this decoupling between SST and q at the mesoscale.

Finally, the last four green markers in Fig. 4 represent s_t in the different domains. As expected, they all are positive: higher SST implies higher surface air temperature. In addition, s_t is strongest in Tradewind (∼0.25) where the SST-associated wind speed and specific humidity variations are weaker. This might result from a decreased advection of other air masses with different temperature which could locally modify temperature values close to the surface.

4.3 LHF Sensitivity to SST and Currents

To evaluate the representation of LHF sensitivity to SST in the EURECA simulation, we perform a linear regression between LHF variations and the SST correction (ΔSST, see Eq. 3) for the different LHF subsets detailed in Sect. 3.2. The results are presented in Fig. 5 for the EURECA domain (first row), the Amazon sub-region (second row), the Downstream sub-region (third row), and the Tradewind sub-region (fourth row).

Figure 5 LHF sensitivity to ΔSST (see Eq. 3) across different regions. Each row corresponds to a specific domain: the EURECA region (first row), the Amazon sub-region (second row), the Downstream sub-region (third row), and the Tradewind sub-region (fourth row). Panels (a), (c), (e), and (g) display the difference between LHF_HR and LHF_LR as a function of ΔSST. Panels (b), (d), (f), and (h) show the differences between LHF_therm and LHF_LR, and between LHF_therm−U and LHF_LR, as a function of ΔSST (blue and orange, respectively). The methodology used to compute these LHF datasets is detailed in the main text. The SST anomaly values are binned into 2 % percentile intervals. Vertical error bars indicate the standard deviation from the mean in each interval, while straight lines represent the least-squares regression fits. The regression slope ± standard error (p value) is reported at the bottom of each panel.

The estimated LHF change per °C associated with the presence of the mesoscale ocean is shown in orange in Fig. 5a, c, e, and g, representing the difference between LHF_HR and LHF_LR. Recall that LHF_HR is computed using the statistically downscaled WRF variables (Eq. 2) in COARE3.5, while LHF_LR corresponds to LHF derived from the smoothed variables, obtained by applying a Gaussian low-pass filter (250 km cutoff length) to the raw data.

The mean LHF sensitivity to mesoscale SST in the EURECA region is 47.7 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ and the Amazon and Downstream sub-regions show similar values (46.7 and 46.9 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ respectively). On the contrary, the sensitivity is lower in the Tradewind subdomain (35 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$). Given the mean LHF values of 150 W m⁻² in EURECA, 130 W m⁻² in Amazon, 160 W m⁻² in Downstream, and 180 W m⁻² in Tradewind, the slopes of the linear regressions correspond to 31.8 % K⁻¹ ($\mathrm{47.7}/\mathrm{150}$), 35.9 % K⁻¹ ($\mathrm{46.7}/\mathrm{130}$), 29.3 % K⁻¹ ($\mathrm{46.9}/\mathrm{160}$), and 19.4 % K⁻¹ ($\mathrm{35}/\mathrm{180}$), respectively. These values align with the theoretical estimate of approximately 33 % K⁻¹ of Fernández et al. (2023). In addition, they are consistent with satellite (33.8 % K⁻¹) and reanalysis (26.6 % K⁻¹) estimates obtained in Fernández et al. (2023) for the whole EURECA domain. The observed regional differences were also reported in that study, with stronger LHF variations near the South American coast (Amazon and Downstream) compared to the open ocean (Tradewind).

We compare the LHF sensitivity results with in-situ observations as well. In Fernández et al. (2024), the LHF changes were assessed across SST mesoscale anomalies of 2 and −0.4 °C in the Amazon and Downstream sub-regions, respectively. The first SST anomaly induced a LHF difference of approximately 160 W m⁻² between itself and its environment (∼80 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$). The second SST anomaly, resulted in a LHF difference of 95 W m⁻² (∼38 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$). Therefore, the latter agrees better with the model estimate for the Downstream subdomain, while the former is significantly larger, likely due to its proximity to the coast. Indeed, the strongest SST gradients in the model also appear over the continental shelf (Fig. 3a) and are not included in this analysis. We should also keep in mind that LHF sensitivity to SST across fronts in in-situ observations is subject to several uncertainty sources which could modify its value and make the direct comparison with model estimates less straightforward. The relative orientation between the sampling device's trajectory and the SST front or even just defining the front's location are among these uncertainty sources.

To quantify the thermodynamic contribution to LHF sensitivity (i.e. the component linked only to SST changes, isolated from the difference between LHF_therm and LHF_LR), we compute its linear regression shown in dark blue in Fig. 5b, d, f, and h. In the EURECA domain, this contribution accounts for a LHF change of 6.6 % K⁻¹ ($\mathrm{9.9}/\mathrm{150}$). In the other sub-regions, it represents 6 % K⁻¹ (Amazon, $\mathrm{7.76}/\mathrm{130}$), 5.2 % K⁻¹ (Downstream, $\mathrm{8.4}/\mathrm{160}$) and 5.3 % K⁻¹ (Tradewind, $\mathrm{9.5}/\mathrm{180}$). Note that in Tradewind, the thermodynamic contribution relative to the total LHF sensitivity to SST is larger than in the other two subdomains. The increasing importance of the thermodynamic contribution towards the open ocean is consistent with the particular air-sea coupling in this region (lower s_u and larger s_t and larger s_q, see Fig. 4). Overall, these LHF sensitivity values show that in the presence of mesoscale anomalies, LHF variations at the mesoscale are mainly associated with the dynamic contribution (i.e. LHF variations linked to the SST-induced modification of the near-surface winds and relative humidity) whereas the thermodynamic contribution remains a second-order contributor, in agreement with previous observation and reanalysis-based studies (Fernández et al., 2023).

For completeness, we also compute ${\text{LHF}}_{\text{therm}-U}-{\text{LHF}}_{\text{LR}}$ as a function of ΔSST, shown in orange in Fig. 5b, d, f, and h. This allows us to separate the surface wind contribution from the relative humidity's within the dynamic component of LHF sensitivity (the air temperature contribution is negligible when compared to the others, not shown). In EURECA, Amazon and Downstream the effect of surface wind variations adds around 8 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ to the thermodynamic contribution. Regarding Tradewind, the wind speed contribution is more modest, representing only an additional 4.7 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$. Therefore, within the dynamic contribution, the majority of the LHF change is linked to relative humidity variations. Indeed, due to the weak SST-q coupling at the mesoscale, q does not follow the Clausius–Clapeyron relation and thus differs from the saturation specific humidity triggering large relative humidity variations and the LHF changes described here. The air temperature effect is even smaller than the surface wind's and it is not shown here.

To quantify the impact of using relative winds (i.e. surface winds minus currents) in LHF computations, we perform a linear regression between LHF_orig−LHF_U and the difference between the norms of relative winds and surface winds. Here, LHF_U represents the dataset computed using only the winds from the first model vertical level (without currents) and LHF_orig stands for the LHF computed using relative winds. The regressions are shown in dark blue in Fig. 6 for the EURECA (panel a), Amazon (panel b), Downstream (panel c), and Tradewind (panel d) regions.

Figure 6 LHF sensitivity to surface currents for (a) EURECA, (b) Amazon, (c) Downstream, and (d) Tradewind. In all panels, blue markers represent the difference between LHF_U and LHF_orig as a function of the difference between the norm of the first vertical level wind velocity and the relative wind velocity field. Orange markers indicate the difference between LHF_orig and LHF_no-CFB as a function of the difference between the norms of relative winds, with and without the CFB effect. The methodology used to derive these LHF datasets is detailed in the main text. The x axis variable is divided into intervals containing an equal number of values, using a 2 % percentile separation. Vertical error bars indicate the standard deviation relative to the mean for each interval, while straight lines denote the least-squares regression fits. The regression slope ± standard error (p value) is displayed in each panel's legend.

Across the EURECA domain, surface wind speed generally exceeds relative wind speed (see Fig. D1a of Appendix D). Winds predominantly blow towards the southwest (−110° from north), while surface currents are oriented northwestward (−70° from north). As a result, considering only surface winds increases LHF by up to 10 W m⁻², a trend observed across all sub-regions. The sign and magnitude of LHF_orig−LHF_U depend on the relative orientation of winds and currents. In the Downstream sub-region (Fig. 6c), the southern edge of an anticyclonic eddy (Fig. 3d) aligns surface currents with surface winds (−120° from the north), leading to $\left|\mathit{U}\right|>|\mathit{U}-{\mathit{U}}_{\mathbf{o}}|$ and thus reducing LHF when relative winds are used. Note that in this region surface currents are the strongest, reaching mean values of 0.45 m s⁻¹ in mid-February (Fig. D1c). On the contrary, Amazon and Tradewind (Fig. 6b and d) exhibit both positive and negative LHF differences, since surface currents and winds are aligned or not depending on the location within the region (Fig. 3c and d). In Amazon, surface current variability is stronger than in Tradewind and covaries with relative wind variations, especially by the end of February (Fig. D1b). On the contrary, the Tradewind subdomain exhibits an enhanced wind variability (values ranging from 4 to 12 m s⁻¹, see Fig. D1d), which mostly drives relative wind variations. In addition, Tradewind surface currents are two times smaller than Amazon's or Downstream's. These LHF variations align well with observations. Using the same representative LHF values for each sub-region, we find that the regression slopes correspond to 10.9, 11.3, 9.4, and 13 % K⁻¹ in the EURECA, Amazon, Downstream, and Tradewind regions, respectively. In-situ observations report a sensitivity range of 5 % K⁻¹–15 % K⁻¹ (Fernández et al., 2024).

Finally, to isolate the CFB contribution to LHF variations, we perform a linear regression between LHF_orig−LHF_no-CFB and the difference between relative wind norms with and without current feedback. The results, shown in orange in Fig. 6, exhibit a similar sensitivity than the regression concerning LHF_orig−LHF_U. Whereas $|\mathit{U}-{\mathit{U}}_{\mathrm{o}}|$ is typically smaller than $\left|\mathit{U}\right|$, the relative velocity in presence of CFB is typically larger than the relative velocity without CFB. In addition, the CFB effect is one order-of-magnitude smaller (within ±3 W m⁻²), consistent with in-situ results from Fernández et al. (2024). According to CFB, surface currents generate wind anomalies in their own direction. Thus, in the Downstream region, where surface currents mostly align with wind direction, CFB enhances relative winds and LHF. In regions without a dominant current direction, LHF variations can be either positive or negative.

4.4 The Amazon Sub-region Vertical Atmosphere and Ocean Structures and Mixed Layer Heat Budget

To delve into the mechanisms behind the modeled mesoscale LHF variations, we examine the atmospheric response configuration associated with SST mesoscale anomalies. Note that the Amazon sub-region is the only of the three considered in this study crossed by the Amazon plume (SSS<35 psu, Fig. 3d). Therefore we focus on it, since we also aim to link the SST mesoscale anomalies to the presence of the Amazon plume and its heat budget. Results for the Downstream and Tradewind sub-regions and EURECA yield similar conclusions in terms of physical mechanisms at play, and, for brevity, we do not include them in this study. The atmospheric analysis follows what is presented in Borgnino et al. (2025) for an atmosphere-only model forced with realistic SSTs.

4.4.1 Atmosphere Vertical Structure

Figure 7 presents the binned distribution of the marine atmospheric boundary layer (MABL) vertical structure, as well as the dependence of SST and SSS on the SST mesoscale anomaly. Over cold SST anomalies, the MABL exhibits increased stability, characterized by positive values of the Brunt–Väisälä frequency anomalies ($N^{\mathrm{2}}=\frac{g}{\mathit{\theta }}\frac{\partial \mathit{\theta }}{\partial z}$, with θ potential temperature), particularly above the MABL height (black line in Fig. 7a). Conversely, over warm SST anomalies, the atmosphere becomes more unstable, with negative values of N².

Figure 7 Panels (a–f) show the vertical structure of the first 2000 m of the atmosphere (in terms of mesoscale anomalies) as a function of the SST mesoscale anomaly: (a) Brunt–Väisälä frequency, (b) horizontal wind speed, (c) saturation specific humidity, (d) specific humidity, (e) specific humidity deficit, and (f) cloud water mixing ratio (QCLOUD). In all cases, the stippling indicates bins whose mean is not significantly different from zero (bin standard deviation larger than the absolute value of the bin mean). In all these panels, the black line represents the MABL height (MABLH), while the purple line indicates LHF. The standard deviation around the MABLH mean is marked in dark shading whereas the standard deviation range for LHF has been omitted since it is much smaller than LHF variations. Panel (g) displays SST (red) and panel (h) SSS (blue) both as functions of the SST mesoscale anomaly. In both cases, the shading represents the ± standard deviation around the bin mean. All calculations are based on daily averages of the model output for JF 2020. Bins are defined using a 2 % percentile range. The values shown here refer to the Amazon sub-region.

Changes in atmospheric vertical stratification also influence wind speed variations (Fig. 7b). Two distinct anomaly dipoles are observed at the warmest and coldest ends of the histogram. Over the highest SST anomalies, wind speed increases slightly near the surface (around 0.05 m s⁻¹), while negative wind speed anomalies dominate above. This behavior aligns with the downward momentum mixing (DMM) mechanism: downward fluxes of momentum from the free troposphere reduce wind speed at higher levels while enhancing it near the surface. The opposite pattern is observed over the coldest SST anomalies in the Amazon sub-region, where momentum transfer towards the surface does not occur: wind speed decreases near the surface and increases at and above the MABL height. Note that PA might as well produce this dipolar pattern in the wind speed histogram. Indeed, Small et al. (2003) showed that in the equatorial Pacific pressure extrema associated with the SST anomalies linked to tropical instability waves formed downwind of the SST extrema as they were advected by the mean flow. While this configuration provides surface wind maxima (minima) over the warmest (coolest) SST mesoscale anomalies, it also results in the formation of moisture and air temperature maxima/minima downwind of the SST maxima/minima. As shown in the following paragraphs, this is not the case here, since air temperature maxima/minima are collocated with SST mesoscale anomaly maxima/minima and air specific humidity does not significantly vary along SST mesoscale anomalies.

Furthermore, the warming (cooling) induced by SST mesoscale anomalies extends well above the MABL height (MABLH) (Fig. 7c). This is reflected in the two-dimensional saturation specific humidity anomaly histogram, where values exceed ±0.1 g kg⁻¹ in the warmest (coldest) intervals. In contrast, specific humidity variations are weaker. Consistent with the Amazon s_q shown in Fig. 4, q slightly decreases over the warmest SST mesoscale anomalies and slightly increases in the MABL located over the coldest SST anomalies. At the MABLH level, we find slightly positive q anomalies which then become negative further above, except for the warmest SST mesoscale anomalies. This pattern of specific humidity anomalies aligns again with the DMM: warmer SSTs destabilize the MABL leading to the entrainment of drier (and colder) air from the free troposphere towards the surface, decreasing surface q (and air temperature, thus reducing s_t with respect to Tradewind where s_q is positive as shown in Fig. 4). Consequently, the specific humidity deficit (Δq) anomalies (Fig. 7e) mostly resemble those of the saturation specific humidity (q_s) in Fig. 7c, with negative Δq anomalies coinciding with negative SST anomalies, and vice versa. The only exception corresponds to the warmest SST mesoscale anomalies where negative anomalies of Δq emerge over the MABLH associated with the reduced q.

It is worth recalling that the weak variations with SST of q at the mesoscale trigger moisture undersaturation imbalances (or relative humidity changes) which are the main driver of mesoscale LHF changes when considering its sensitivity to SST (Fig. 5). If DMM were not present, there would not be a downward momentum flux of drier air from the top of the MABL that reduces the specific humidity variations close to the surface. We would then expect higher (lower) values of q over the warmer (cooler) SSTs as predicted by the Clausius–Clapeyron equation, thus reducing the relative humidity variation range. This is the situation found for scales larger than the mesoscale, where air over the ocean is almost saturated. Such a configuration would reduce the specific humidity deficit modulations shown in Fig. 7e and therefore LHF changes.

Figure 7f presents the cloud water mixing ratio (QCLOUD). Although anomalies remain weak, an increase in QCLOUD is observed at the MABLH level for the majority of the SST mesoscale anomaly range. In regions of negative SST mesoscale anomalies, the strongest QCLOUD increases are confined close to the MABLH. This finding aligns with previous observational (Acquistapace et al., 2022) and modeling (Borgnino et al., 2025) studies, which detected stratiform shallow clouds over stable tropical MABLs. Over positive SST mesoscale anomalies, QCLOUD positive anomalies extend upward potentially yielding more vertically developed clouds in the context of a more unstable atmosphere.

In all the panels discussed above, the purple line represents the mean binned LHF as a function of the SST mesoscale anomaly. Consistent with observational results from both satellite (Fernández et al., 2023) and in-situ data (Fernández et al., 2024), a significant LHF increase of approximately 50 W m⁻² is observed. This increase is driven by the SST-induced modifications of saturation specific humidity maintaining the smoothed relative humidity (thermodynamic contribution, quantified with LHF_therm−LHF_LR) and the modulations of relative humidity and wind speed (dynamic contribution, the contribution of wind speed to the dynamic contribution is quantified with ${\text{LHF}}_{\text{therm}-U}-{\text{LHF}}_{\text{LR}}$) induced by the SST mesoscale anomalies. All of them were assessed separately in the previous section. Additionally, the black line represents the simulated MABLH mean per SST mesoscale anomaly bin. Its values remain around 600 m for most SST anomaly bins, consistent with expected tropical ocean conditions. Thus, the MABL response does not fully align with the downward momentum mixing (DMM) mechanism, although changes in vertical stratification, wind speed and specific humidity suggest that DMM is active at the mesoscale in the Amazon sub-region. We hypothesize that this discrepancy may arise from the MABL stability parameterizations used within the model. A more detailed analysis would be required to confirm this hypothesis, but it is beyond the scope of this study.

Note that all the anomalies (atmospheric and SST) shown in Fig. 7 are small compared to the ones shown in Borgnino et al. (2025). A direct comparison is not straightforward since their simulation is SST-forced and ours is coupled, and they use a different filtering procedure: no time filter and a high-pass Gaussian filter with a 150 km cutoff length. In addition, they evaluate the vertical atmospheric structure of the whole EURECA domain (not only Amazon). However, in the range of the anomalies presented in this article (±0.3 °C), the results are consistent. For example, in their Fig. 4d wind speed mesoscale anomalies range from −0.2 to 0.2 m s⁻¹ and SST mesoscale anomalies from −0.6 to 0.5 °C. Nevertheless, between ±0.3 °C their wind speed anomaly values lie within the ±0.05 m s⁻¹ range of Fig. 7b.

Finally, Fig. 7g and h illustrate the dependence of SST and SSS on SST mesoscale anomalies. As expected, SST increases monotonically with the SST mesoscale anomaly: the warm SST anomalies are associated with higher SST values. This is not the case for SSS. It shows a minimum around a SST mesoscale anomaly of 0.1 °C. However, the spread of SSS values along bins is large. This means that a given SST mesoscale in the binned distribution might contain SSS data from different locations within the Amazon sub-region, which prevents us from directly associating the SSS minimum to the exclusive presence of the Amazon plume.

4.4.2 Ocean Surface Features

To gain a deeper insight into the relation between the SST and the Amazon plume, Fig. 8 presents four snapshots of the SST and surface current fields in February 2020 in the Amazon sub-region. In the beginning of February 2020 (Fig. 8a), the southern end of the cold filament (Fig. 3a), characterized by SST<25.8 °C, is advected northward and westward by the surface circulation. It is not until the 16 February that the Amazon plume, delineated by the 35 psu isoline (Reverdin et al., 2021), reaches the southeastern end of Amazon (Fig. 8b). For this reason, the remaining analyses concerning the presence of the Amazon plume focus only on February 2020 instead of the whole JF season. At this stage, the Amazon plume is mostly warmer than its environment (>26.9 °C) although it starts mixing with the cold filament on its southwestern flank. This results in SST values below 26 °C within the Amazon plume boundaries. Finally, some of the plume-related warmest waters are advected northward out of the 35 psu isoline.

Figure 8 Daily snapshots of SST (shading) overlaid with surface currents (arrows) and the 35 psu isoline which we use to delimit the Amazon plume (magenta) for the (a) 9, (b) 16, (c) 21, and (d) 27 February 2020.

By the 21 February 2020 (Fig. 8c), the Amazon plume arrives in the northern part of Amazon and preserves warmer waters than its environment, still reaching 26.9 °C. However, colder waters resulting from the horizontal mixing with the cold filament still remain inside the limits of the plume on its southwestern part. Finally, by the end of the month (Fig. 8d), the Amazon plume is advected westward and northward, following the rotation of a cyclonic eddy whose center lies at 8° N, 55.25° W. It loses its shape and the SST field becomes more homogeneous within the whole Amazon sub-region.

Therefore, during several days in February 2020, a large fraction of the Amazon plume surface presents warmer waters than its environment (>26.9 °C). These SST values correspond to positive SST mesoscale anomalies as shown in Fig. 7g. Hence, the Amazon plume acts on LHF as any other non-plume warm SST mesoscale anomaly would (Figs. 5 and 7), via both the thermodynamic and dynamic contributions assessed in previous sections. The same applies in the opposite direction for the colder southwestern part of the plume which interacts with the cooler coastal band, whose SST (<26 °C) would belong to the negative SST mesoscale anomaly bins of Fig. 7g.

We explored other plume-related mechanisms which affect LHF, such as the potential subsurface heat release detected in in-situ observations (Fernández et al., 2024). If the base of the isothermal layer (THERM) is deeper than the mixed layer depth (MLD), subsurface warm inversions can develop and potentially be released to the atmosphere mainly in the form of LHF provided that salinity-driven stratification is reduced. We checked this and found that, contrary to in-situ observations, the model does not represent such inversions. Consequently, their contribution to LHF spatial variability cannot be isolated.

4.4.3 Vertical Ocean Structure

The two parts of the Amazon plume in terms of SST are also observed when averaging in February 2020: a warmer core on its eastern part and a cooler area on its southwest linked to the coastal filament (Fig. 9a). The warmer SST part belongs to a larger southeastern-northwestern warm SST band which is advected through the Amazon subdomain from the southeast to the northwest by surface currents. It coincides with a sharp SSS gradient which marks the transition between the fresh plume waters to the saltier open-ocean side of Amazon. This band is characterized by shallow MLDs (less than 15 m). Cooler and saltier waters dominate towards the northeast of Amazon and to its northeast, together with thicker MLs. Note that the MLDs are larger in the east (>50 m) than in the west (<40 m) even if they have similar SST (∼26.4 °C) since salinity is higher in the east (>36 psu). Note that colder southwestern area is characterized by shallower MLDs than the northeastern most open-ocean part of Amazon despite being colder. This is due to the mixing with the Amazon plume which reduces salinity resulting in thinner MLs.

Figure 9 February 2020 average of (a) SST (shading) and surface currents (arrows), (b) mixed layer depth (MLD, shading), (c) barrier layer thickness (BLT, shading), and (d) OSS index integrated down to the base of the mixed layer (〈OSS〉, shading). In all panels black contours denote the February 2020 averaged SSS and the magenta contour the mean February 2020 position of the edge of the Amazon plume (35 psu isoline).

Additionally, we present the spatial distribution of the barrier layer thickness (BLT) and the OSS index integrated down to the ML base with superimposed SSS contours in Fig. 9c and d, respectively. The BLT map clearly shows that to the west and northwest of the plume, the mixing of cooler pre-plume waters with the warm and fresh waters preceding the Amazon plume triggers a decoupling between haline and thermal stratifications and leads to the formation of barrier layers (BL) as thick as 35 m. In the rest of the domain including the interior of the Amazon plume, the MLD and the BLT are close to each other. Although this configuration mostly appears in the most open-ocean region of Amazon, where waters are the saltiest, it is also present within the Amazon plume and prevents the formation of subsurface temperature inversions there. Finally, the OSS index map shows that, despite the SSS heterogeneity of the region, temperature dominates over salinity in driving total stratification everywhere (〈OSS〉<50 %). This configuration is more pronounced in the most open-ocean part of Amazon. Only some regions to the northwest of the Amazon plume and in the southwest of the Amazon plume show 〈OSS〉 values higher than 30 %. The influence of the Amazon plume being advected toward these areas before temperature can adjust might explain these values: salinity rapidly decreases whereas temperature remains temporarily constant.

4.4.4 Mixed Layer Heat Budget

To further investigate heat transfer in and out of the ML, Fig. 10 displays the various terms of the ML heat budget (Eq. 4) averaged over February 2020. The total temperature tendency (∂_tT) exhibits a highly heterogeneous pattern (Fig. 10a). In the southwestern part of the domain, warming is primarily driven by the horizontal advection of warmer waters from the east (Figs. 10b and 9a). Some of these waters lie within the boundaries of the Amazon plume. The warming effect of horizontal advection is partly compensated by the cooling contribution of the residual (Fig. 10f), which can be mainly attributed to vertical diffusion in this region as it meets the Cronin et al. (2015) criteria (detailed in Appendix B).

Figure 10 Averaged February 2020 mixed layer heat budget (terms of Eq. 4) vertically integrated down to the mixed layer depth for the Amazon sub-region. (a) Total temperature tendency, (b) horizontal advection, (c) vertical advection, (d) entrainment, (e) atmospheric forcing, and (f) residual. The hatched areas in (f) indicate locations where the residual is primarily associated with vertical diffusion, following Cronin et al. (2015). In all panels, black contours denote averaged February 2020 SSS and the magenta contour represents the mean February 2020 position of the edge of the Amazon plume (35 psu isoline).

In contrast, the Amazon plume-related waters located in the east and northern parts of the plume experience net cooling, mainly due to atmospheric forcing (Fig. 10e). An evaluation of the magnitude of the different radiative and turbulent air-sea heat fluxes contained in the atmospheric forcing term (Eq. 4) shows that the enhanced LHF over the warmer waters of this region is the main driver of this heat loss (not shown). The thinner MLs in this region also facilitate the enhanced heat transfer to the atmosphere (Fig. 9b). In addition, the nearly uniform temperature in the warmest eastern part of the plume (south of 7.75° N and east of 54.25° W) explains the absence of strong horizontal advection values in its interior (Fig. 10b), as warm water from the surroundings is advected into a region where little temperature variation occurs. Finally, the warmer waters (SST>26.7 °C) located over the sharp salinity gradient (SSS between 35–35.7 psu) exhibit a slightly positive ∂_tT. This increase results primarily from horizontal advection from the southeast which transports warm plume-related waters northwestward.

For completeness, Fig. 10c shows the vertical advection contribution to the total temperature tendency. It is one order of magnitude smaller than the horizontal advection everywhere, not significantly contributing to the total temperature tendency. The same applies to the entrainment term (Fig. 10d), which depicts slightly positive entrainment values in the plume's eastward-northeastward flank, likely resulting from MLD variations associated with the lateral displacement of the plume. Another contributing factor to ∂_tT is vertical diffusion. Positive values of the mixed layer heat budget residual appear around the 35.5 to 35.9 psu isolines, to the east and northeast of the plume (Fig. 10f). This residual can primarily be attributed to vertical diffusion, as it meets the Cronin et al. (2015) criteria.

All the processes described above result in a weakening of the SST anomalies associated with the Amazon plume and therefore, an homogenization of the SST field in the Amazon sub-region as February 2020 draws on. In the southwestern part of the plume, where the interaction with the cold coastal filament occurs and SST is the lowest, temperature in the mixed layer increases due to warm horizontal advection from the east, where SST presents its maximum. On the contrary, in the center and eastern parts of the plume, where SST is the highest, there is a cooling tendency mainly associated to heat loss to the atmosphere, mostly in the form of LHF (not shown).

5 Discussion and Conclusion

High-resolution coupled simulations serve as a powerful complement to remote sensing, reanalysis, and in-situ observations to study ocean mesoscale influences on latent heat flux (LHF) variability in the Northwest Tropical Atlantic. In this study, we employ the WRF-CROCO coupled EURECA simulation at 1 km oceanic and 2 km atmospheric resolution, fully resolving ocean mesoscale (O(50–250 km)) processes. We focus on January–February 2020, when the Intertropical Convergence Zone (ITCZ) shifts southward, and analyze four domains: the full EURECA region; Amazon and Downstream (coastal regions, with the former influenced by the Amazon plume); and Tradewind (a more quiescent open-ocean region).

Our analysis of air-sea coupling coefficients confirms a robust mesoscale positive correlation between surface current vorticity and wind curl across all regions. The associated coupling coefficient is denoted as s_w. This correlation is the strongest in the Downstream sub-region (s_w=0.26, Fig. 4e) compared to Tradewind (s_w=0.24) and Amazon (s_w=0.21). Although these variations might seem small, the error bars associated with these coupling coefficients do not overlap (Fig. 4, red markers) and all s_w values are statistically significant (see Appendix C). They represent changes of 14.3 % in s_w: (0.26–0.22)$/\mathrm{0.28}$. This difference might be related to the presence of stronger surface currents in Downstream than in Amazon and Tradewind associated with the nearly-stationary eddy located in front of Trinidad and Tobago (Fig. 3d). Stronger currents may drive a more robust wind response, less likely to be masked by other atmospheric processes.

Most coupling coefficients at the mesoscale are weaker than satellite-derived values in Fernández et al. (2023) but agree with those reported in Renault et al. (2019b). We find s_u=0.26 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ for the EURECA domain, peaking in Amazon and Downstream (0.29 and 0.31 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ respectively, Fig. 4) and reaching its minimum in Tradewind (0.11 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, Fig. 4). This pattern reflects stronger Downward Momentum Mixing (DMM) activity in Amazon and Downstream than in Tradewind. Meanwhile, the mesoscale SST-q correlation (s_q) is weaker than the Clausius–Clapeyron estimate (1.3 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) in all cases: −0.05, −0.09, and 0.24 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ in Amazon and Downstream and Tradewind respectively. Finally, s_t is always positive and ranges from 0.16 in Amazon (in the warm eddy corridor) to 0.25 in Tradewind (open ocean).

We quantify LHF sensitivity to sea-surface temperature (SST) anomalies through linear regression analysis (Fig. 5). Consistent with satellite observations (Fernández et al., 2023), LHF increases by 31.8 % per 1 °C warming, with the strongest sensitivity in Amazon (35.9 % K⁻¹, Fig. 5c) and the weakest in Tradewind (19.4 % K⁻¹, Fig. 5g). Furthermore, we separate the thermodynamic (mesoscale SST effects in the saturation specific humidity without accounting for any change in relative humidity) and the dynamic contribution (SST-induced modifications of the near-surface atmosphere manifested in surface wind speed and relative humidity changes), confirming that dynamic effects dominate. The thermodynamic contribution represents only 4.5 % per 1 °C–5.5 % per 1 °C and the rest is related to the dynamic contribution. Within the dynamic contribution, the fact that mesoscale specific humidity does not evolve as predicted by Clausius–Clapeyron produces an undersaturation imbalance (translated into larger relative humidity variations) and drives most of the LHF changes. Finally, we assess LHF sensitivity to surface currents (Fig. 6). Accounting for relative winds instead of absolute winds induces LHF variations up to 15 W m⁻², consistent with in-situ results (Fernández et al., 2024). The current feedback (CFB) effect is much smaller, contributing at most with 3 W m⁻².

The vertical structure of the marine atmospheric boundary layer in the Amazon sub-region reveals key mechanisms behind the observed coupling. We observe DMM-consistent patterns, with a dipole in N² (more unstable over warm SSTs, more stable over cold SSTs) and a reduction (increase) of near-surface specific humidity (near-surface winds) over warm SST anomalies. Similar findings hold for Downstream and Tradewind (not shown). In the ocean, the Amazon plume-related waters are characterized by shallow mixed layer depths (Fig. 9b) and most of them are warmer than their environment (Fig. 8b and c) leading to the formation of warm SST mesoscale anomalies. A smaller fraction of them presents lower SSTs due to the interaction with a cooler coastal filament. Thus, the Amazon plume affects LHF through its associated mesoscale anomalies as described above. No temperature inversions underneath the MLD and potential subsurface warm layer heat release in the borders of the plume was detected in this simulation.

Finally, we perform a mixed layer heat budget analysis in the Amazon sub-region (Fig. 10) and find that the total temperature tendency (∂_tT) acts to homogenize the SST field between the inside and outside of the plume as time draws on. In particular, the warmest part of the plume experiences net cooling due to enhanced heat loss to the atmosphere whereas the coolest region of the plume becomes warmer as a consequence of horizontal advection from the warmest part of the plume.

This study extends the mesoscale findings of Fernández et al. (2023, 2024) providing a regionalized analysis of the air-sea coupling and the ocean-atmosphere vertical structure. It also links the air-sea coupling with the Amazon plume. However, our analysis is restricted to boreal winter (January and February). Future work should investigate other seasons in the EURECA simulation (June 2019–June 2020), as coupling coefficients vary seasonally (Conejero et al., 2024), with stronger air-sea interactions in boreal winter than in summer. Additionally, the Amazon runoff peaks in summer, amplifying the Amazon plume's influence on LHF heterogeneity. Moreover, this study does not address the ocean submesoscale (O(<50 km)) although evidence suggests that submesoscale ocean structures also impact the near-surface atmosphere (Meroni et al., 2018; Gaube et al., 2019). This is mainly because daily averaging, used to remove the diurnal cycle, might also eliminate important submesoscale signals. A more sophisticated approach, such as multichannel singular spectrum analysis (M-SSA), could better preserve submesoscale variability. Furthermore, comparing the higher-resolution EURECA simulation (1 km ocean, 2 km atmosphere) with the coarser Antilles simulation (2.5 km ocean, 6 km atmosphere) would clarify the impact of spatial resolution on mesoscale and submesoscale air-sea interactions. Finally, longer simulations (≥10 years) would enable interannual variability studies, particularly regarding Amazon plume detachment, which occurs irregularly (i.e. absent in 2010, 2011, and 2013; Olivier et al., 2022).

Our findings emphasize the need for high-resolution modeling in climate studies. Traditional climate models rely on coarse SST grids, which suppress the small-scale air-sea disequilibrium that governs LHF release. Implementing the LHF downscaling algorithm used here, which has been shown to improve LHF representation by a factor of two in SST-forced WRF simulations (see Fig. 8b of Fernández et al., 2023) in model couplers could enhance air-sea flux estimation in high-resolution climate simulations. Future work should explore its implementation within global coupled models to improve energy exchange parameterizations.

Appendix A: Numerical Considerations in the Mixed Layer Heat Budget Computation

In the EURECA simulation, only the monthly means of the advective-diffusive equation terms are stored. For temperature, this equation reads:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }\underset{\text{Total tendency}}{\underbrace{{\partial }_{t}T}}+\underset{\text{Advection}}{\underbrace{\mathrm{\nabla }\cdot \left(\mathit{u}T\right)}}=\underset{\text{Forcing}}{\underbrace{\mathcal{F}}}+\underset{\text{Mixing}}{\underbrace{\mathcal{D}}}.\end{array}$$

One of the challenges in computing the mixed layer (ML) heat budget is ensuring its closure, as numerical choices in the calculation of each term can introduce discrepancies. To mitigate this issue, we followed the computational steps detailed in the CROCO documentation when computing all the terms in Eq. (A1) on a daily basis.

Figure A1 presents a comparison between the recomputed and stored Eq. (A1) terms, integrated down to 100 m (a depth beyond any mixed layer depth to minimize entrainment effects). The left column (except for the last row) displays the recomputed terms, the central column shows the corresponding stored values from the simulation output, and the right column represents the numerical bias, i.e. the difference between the recomputed and stored values.

Figure A1 Comparison between the stored advective-diffusive equation terms (Eq. A1) and the recomputed ones, integrated over a layer of 100 m, (deeper than the maximum MLD recorded).

The total temperature tendency, averaged over February 2020, exhibits a consistent spatial pattern between the recomputed (Fig. A1a) and stored (Fig. A1b) values. A region of strong cooling extends from the southern part of the domain (around 54° W) towards the north and northwest, surrounded by positive temperature tendency values, particularly in the southernmost part of the domain. The numerical bias (Fig. A1c) is one to two orders of magnitude smaller than the actual temperature tendency values, indicating that the numerical error in recomputing the total tendency is relatively small.

However, significant differences arise when comparing the advection terms. The recomputed horizontal (Fig. A1d) and vertical (Fig. A1g) advections differ substantially in both pattern and intensity from their stored counterparts (Fig. A1e and h, respectively). The recomputed advections are systematically stronger, and their differences (Fig. A1f and i) are of the same order of magnitude as, or even exceed, the values themselves.

In contrast, the atmospheric forcing term remains consistent between the recomputed (Fig. A1j) and stored (Fig. A1k) datasets, with negligible numerical bias (Fig. A1l). This suggests that numerical errors as well as high frequencies (timescales faster than 1 d) primarily impact the advection terms, while the atmospheric forcing is well represented in the recomputed dataset.

The discrepancy in the advection terms leads to a significantly large residual field (Fig. A1m), which exceeds the total temperature tendency itself (Fig. A1a). This stands in contrast to the stored residual field (Fig. A1n), where significant values are only observed around 7.5° N. The difference between the recomputed and stored residual fields (Fig. A1o) is almost as large as the recomputed residual field itself, highlighting the impact of numerical discrepancies. Additionally, the stored residual field is evenly partitioned between horizontal and vertical mixing, whereas these two contributions cannot be separately estimated from the recomputed fields.

Overall, Fig. A1 underscores the critical role of numerical schemes in computing heat budgets, particularly for advection and derivative calculations. In this case, the recomputed horizontal and vertical advections were obtained using a second-order centered scheme, whereas the simulation employs a fifth-order upstream advection scheme, leading to significant differences in the resulting fields. Thus, in the mixed-layer heat budget-related results, we will work with the stored advection terms whereas the rest of the terms (total tendency, entrainment and atmospheric forcing) will be recomputed each day from model data and then averaged in a monthly basis.

Appendix B: Criteria to Associate the Residual from the Mixed Layer Heat Budget to Vertical Diffusion

A common approach to estimating vertical heat diffusivity is to infer it from the residual of the mixed layer heat budget (Girishkumar et al., 2020). Cronin et al. (2015) extended this approach by proposing selection criteria for estimating vertical heat diffusivity at the ML base. First, if the vertical temperature gradient is small, the residual may not be associated with vertical heat diffusivity. Therefore, we require a vertical temperature gradient greater than 0.0003 °C m⁻¹ within 5 m below the MLD for the residual to be considered representative of vertical diffusion. However, as noted by Cronin et al. (2015), in the presence of strong currents, the residual term is primarily influenced by heat flux convergence in a stratified shear flow. To minimize this effect, we impose an additional filtering criterion where the residual is considered representative of vertical diffusion only when horizontal advection remains within two standard deviations of its mean value.

The formation of a barrier layer (BL) inhibits heat mixing below the ML. When the BL is thin (less than approximately 15–20 m), the MLD and the isothermal layer (THERM) are close to one another, typically resulting in a negative temperature gradient at the ML base. In contrast, when the BL is thick, the base of THERM is significantly deeper than the MLD, leading to a more uniform temperature profile within the BL, which can result in small temperature gradients or even temperature inversions. The barrier layer thickness (BLT) is defined as the difference between THERM and MLD, when THERM is deeper than MLD (Sprintall and Tomczak, 1992). Here, we define THERM as the deepest level at which temperature decreases by at least 0.2 °C relative to the 10 m-depth temperature, following Gévaudan et al. (2021). To ensure that the residual reflects vertical diffusion, we impose a final criterion: the BL must either be thinner than 15 m or contain a temperature inversion of magnitude greater than 0.2 °C.

In summary:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{Residual represents vertical diffusion where}\\ \text{(B1)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}T}{\mathrm{d}z}>\mathrm{0.0003}\mathrm{°}\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}\mathrm{5}\mathrm{m}\text{ below the MLD or}\\ {\mathrm{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{u}}_{\mathrm{h}}\cdot T\right)<\mathrm{2}\cdot {\mathit{\sigma }}_{{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left({\mathit{u}}_{\mathrm{h}}\cdot T\right)}\text{ at the MLD or}\\ \text{BLT}<\mathrm{15}\mathrm{m}\text{ or}\\ \text{BLT}>\mathrm{15}\mathrm{m}\text{ and }\mathrm{\Delta }T>+\mathrm{0.2}\mathrm{°}\mathrm{C}\text{ within the BL}.\end{array}\right.\end{array}$$

In the second condition, the subscript h indicates that only horizontal derivatives are considered.

Appendix C: Coupling Coefficient Linear Regressions

This section presents the binned linear regressions between mesoscale anomalies whose slopes are displayed in Fig. 4: the coupling coefficients. A discussion follows in Sect. 4.2.

Figure C1 Binned scatter plots of surface current vorticity (ξ) vs. surface wind curl mesoscale anomalies in the (a) EURECA, (b) Amazon, (c) Downstream, and (d) Tradewind. In all cases, error bars represent the standard deviation of each bin. Bins are computed using daily averages from the JF 2020 season, excluding the continental shelf (seafloor depth <100 m) and islands as described in the main text. All the panels include least-squares regression lines for the mesoscale, with the slope ± standard error (p value) indicated in the legends. To build each bin, surface current vorticity samples are divided into 2 % percentile intervals.

Figure C2 As in Fig. C1 but for the linear regression between near-surface wind (U) and SST mesoscale anomalies.

Figure C3 As in Figs. C1 and C2 but for the linear regression between near-surface specific humidity (q) and SST mesoscale anomalies.

Figure C4 As in Figs. C1, C2, and C3 but for the linear regression between near-surface air temperature (T) and SST mesoscale anomalies.

Appendix D: Time Series of Surface Winds, Relative Winds and Surface Currents

Figure D1 Time series of the area-weighted mean wind speed (blue), relative wind (red) and surface current (black) for (a) EURECA, (b) Amazon, (c) Downstream, and (d) Tradewind.