The turbulent heat fluxes (THFs) between the ocean and the atmosphere account for the exchange of energy caused by the thermal imbalance between the two fluids (sensible heat flux, SHF), and the exchange of energy due to the sub-saturation and following water phase change at the air-sea interface (latent heat flux, LHF). Both fluxes are generally described with bulk formulations which strongly rely on the Monin-Obukhov similarity theory (MOST) (Monin and Obukhov, 1954). In particular, they are written in terms of mean quantities that include a bulk transfer coefficient, the near-surface wind speed and an air-sea imbalance term: the temperature difference for SHF and the humidity saturation deficit for LHF (Fairall et al., 2003).

An accurate measure and computation of THFs is still a challenging task. For instance, Cronin et al. (2019) stress that the current ocean observing system does not meet the necessary requirements to accurately measure THFs, mainly due to a lack of global coverage of surface humidity and air temperature. Other studies highlight the challenges in representing THFs as bulk quantities using the MOST (Fairall et al., 2003; Edson et al., 2013) as the uncertainties in these parametrizations are large. Moreover, when such bulk formulations are integrated over the ocean basins they cause a large imbalance in the energy budgets (Yu, 2019). As a consequence, there is a far-reaching ongoing effort to properly define and tune these THFs in numerical weather predictions as well as in Earth System models (ESMs). Differences of the order of 15% have been observed between model THF data sets and have been mainly attributed to the choice of parametrization scheme of the bulk formulae. However, other sea surface temperature (SST) and humidity-related approximations, such as the skin temperature correction and the salt-related reduction of the saturation specific humidity, have been found to substantially impact THFs estimates as well (Brodeau et al., 2017).

It is a common practice in the literature to separate between the ocean large-scale and small-scale (mesoscale, O(10-200) km) as their interaction with the atmosphere has been suggested to be different (Chelton and Xie, 2010; Small et al., 2019; Gentemann et al., 2020). At the large-scale, the atmospheric dynamics has been shown to drive ocean variability (Gill, 1982). However, at scales smaller than 200 km, the ocean actively forces the near-surface atmosphere impacting air temperature, frictional stress and the marine atmospheric boundary layer (MABL) stability (Small et al., 2008).

The importance of the ocean dynamics in affecting the SST variability through THFs, the so-called ‘ocean weather’, has been highlighted at the mesoscale, O(10-200) km, on monthly (and longer) time scales (Bishop et al., 2017). Small et al. (2019) exploit observational and numerical simulation data to characterize the monthly and inter-annual THFs variability. They find that both, the atmospheric and the oceanic intrinsic variability, contribute to the heat flux variability, with SST contributing the most in mid-latitude ocean frontal regions whereas the wind influence is prominent in tropical and sub-tropical regions. Mesoscale eddies have been shown to impact surface THFs with effects on the surface wind, cloud cover and rainfall throughout the extra-tropical ocean: in the Gulf Stream region (Minobe et al., 2008), in the Southern Ocean (Frenger et al., 2013), in the Kuroshio extension (Xu et al., 2011; Ma et al., 2015; Chen et al., 2017), in the South China Sea (Liu et al., 2018; Liu et al., 2020), in the Agulhas (O’neill et al., 2005) and Malvinas currents (Villas Bôas et al., 2015; Leyba et al., 2017). Warm eddies are found to induce an increase in surface fluxes, surface winds and vertical motion enhancing cloud cover and precipitation. The dependence of the THF intensity on eddy size has been highlighted by Lin and Wang (2021), with larger eddies driving stronger THFs, and a more intense atmospheric response.

Moreton et al. (2021), using coupled model data, are the first to quantify the turbulent heat flux feedback (THFF) at the ocean mesoscale. The THFF measures the damping rate of SST anomalies as a consequence of the energy gain or loss associated with THFs and its magnitude in models has been found to critically depend on the atmospheric grid spacing. Although these results seem to be highly model-dependent, a better understanding of THFF is key to obtain an accurate estimate of THFs. Strobach et al. (2020) show how the air-sea coupling, in particular in terms of the wind response to the SST forcing and its feedback on the THFs, is responsible for a three-to-six-day oscillation in both surface wind speed and SST.

In the tropical oceans LHF generally increases with SST. However, the connection between THF and the ocean temperatures is subject to many different local and regional factors (Kumar et al., 2017). For instance, the seasonality of the Asian monsoon triggers a larger specific humidity seasonal cycle over the Arabian Sea and the Bay of Bengal when compared to other tropical basins. This is linked to the fact that during winter the monsoonal winds blow dry continental air over the ocean, resulting in large LHF. During summer, the air moving towards the continent is relatively moist when going over the bays, bringing near-surface air close to saturation and reducing LHF despite the increase in SST.

There is evidence that also at the ocean submesoscale, O (1-10) km, SST gradients can affect the surface wind response (Gaube et al., 2019) and can produce THFs that are significantly larger than state-of-the-art bulk parametrization outputs (Shao et al., 2019). Cloud-resolving high resolution numerical simulations also show the importance of submesoscale SST structures in rapidly modulating the surface wind field (Meroni et al., 2018) and in generating strong turbulent air-sea fluxes that drive significant atmospheric responses (Strobach et al., 2022). At larger spatial scales, but still on relatively short temporal scales (days to weeks), Ma et al. (2020) highlight that the cold wakes of tropical cyclones have a similar imprint on the surface fluxes and the lower atmosphere. In fact, the cold wakes reduce the THFs, slowing down the surface wind and reducing cloud fraction and rainfall, by means of a cross-track secondary circulation (Pasquero et al., 2021).

In-situ observations, despite being sparse and intermittent, provide unique and essential data to understand the physical processes at play. The EUREC⁴A (ElUcidating the RolE of Clouds-Circulation Coupling in ClimAte, www.eurec4a.eu) initiative (Bony et al., 2017) aims to advance the understanding of the interplay between clouds, convection and circulation (and their role in climate) with an unprecedented extensive observational effort. The core activity of EUREC⁴A has been a 6-week field campaign between the 12th of January and the 23rd of February 2020 in the north-western tropical Atlantic. High-resolution, synchronized observational data have been collected using cutting-edge technology on airplanes, ships, autonomous vehicles, augmented with the Barbados Cloud Observatory time series (Stevens et al., 2021). In addition to the atmospheric observations, an unprecedented sampling of the upper ocean layers and the air-sea interface has been achieved. This was made possible by additional initiatives such as the European EUREC⁴A-OA (Karstensen et al., 2020; Speich et al., 2021, EUREC⁴A-OA) and the American ATOMIC (Quinn et al., 2021, Atlantic Tradewind Ocean-atmosphere Mesoscale Interaction Campaign) ones. In particular, EUREC⁴A-OA constitutes the ocean component of EUREC⁴A. It investigates heat, momentum, freshwater and CO₂ transport in the ocean (Reverdin et al., 2021), their exchanges across the air-sea interface (Olivier et al., 2022) and the related atmospheric boundary layer processes (Acquistapace et al., 2022). In particular, EUREC⁴A-OA focuses on small-scale ocean dynamics (0.1-100 km). During the campaign, a wide set of innovative and standard observing platforms were deployed, including Saildrones, ocean gliders, wave gliders, different types of surface buoys, profiling floats and 4 research vessels. Motivated by EUREC⁴A-OA, this study focuses on the 5°-17°N, 60°-51°W box, hereby referred to as the EUREC⁴A-OA region, during the December-January-February (DJF) season.

Given the importance of THF variability at the ocean small scales (Shao et al., 2019; Strobach et al., 2022) and in particular the LHF variability, the main goal of the present work is to understand the coupling mechanisms between the small-scale SST features and the near-surface atmosphere as well as their impacts on LHF. The rest of the manuscript is organized as follows. The data employed are described in section 2. The methodology and the theoretical framework linking SST with LHF is derived in section 3. Section 4 contains the main results and is followed by a summary and the main conclusions in section 5.

THFs are often computed using bulk formulae. In particular, a widely-used expression for LHF (Yu, 2009) reads

where ρ_a and L_v represent air density and latent heat of evaporation, respectively. C_e is the moisture exchange coefficient and accounts for the effects of the radiation balance and atmospheric stability on LHF. q_* is the saturation specific humidity, which, over the ocean, is computed using the SST values in the Clausius-Clapeyron equation.

Moreover, this paper also makes use of the COARE3.5 bulk algorithm to obtain LHF from SST, U_{10 m} , 2 m relative humidity (RH _2m ) and T _2m by means of the MOST (Fairall et al., 1996; Weller and Anderson, 1996; Fairall et al., 2003; Edson et al., 2013). COARE3.5 also needs information about the surface radiation balance, surface pressure, reference heights for the large-scale observations, latitude and height of the boundary layer to compute turbulent fluxes. For simplicity, all these parameters are set to the constant reference values presented in Table 1 . No rain rate or wave-related corrections are considered and fluxes are computed from skin temperatures, with the exception of WRF data where bulk SST is considered. Therefore, the cool skin correction to the SST is not applied to SeaFlux and ERA5 and is considered in WRF.

As stated in the Introduction, the relation between atmospheric and near-surface oceanic variables is quite different depending on the spatial scale considered, with the atmosphere generally forcing the ocean at the large scale and the ocean forcing the atmosphere at the mesoscale and below. In order to separate the large from the small-scale patterns, U_{10 m} , q _2m , T _2m and SST are filtered with an isotropic Gaussian filter as described in the Appendix. Hence, a given variable ψ is decomposed as

with ψ ¯ being the filtered large-scale component and ψ ^′ the residual field, which contains small-scale information. The standard deviation of the Gaussian filter σ represents the threshold imposed to separate the small and the large-scale during the filtering (see the Appendix). Unless otherwise stated σ = 150 km, which roughly corresponds to the oceanic mesoscale; it also separates the larger scales at which the atmospheric winds directly impact ocean mixing and SST from the smaller scales at which the upper ocean thermal properties drive an atmospheric response (Seo, 2017; Gentemann et al., 2020).

The residual fields are used to study the strength of the air-ocean coupling. We indeed compute, through a linear regression, the coupling coefficients between SST’ and other variables such as U’_{10 m} , q’_{2 m} and T’_{2 m} . The resulting coupling coefficients are denoted as α_U , α_q and α_T respectively and their statistical significance is assessed using a two-sided t test. These coupling coefficients are computed taking into account the residual fields at all (daily) time steps in a subset of the grid points in the EUREC⁴A-OA region. A sub-sampling has been introduced in order to consider spatially uncorrelated data, as detailed in the Appendix. A generic coupling coefficient is written as

The coupling coefficients are the core of the downscaling linear algorithm proposed in this study since they enclose all air-sea interaction processes. Given a certain SeaFlux/ERA5 variable (ψ) we build a new data set of the same variable including fine-scale SST features (ψ_HR ) as

where ΔSST stands for the SST correction and represents the deviation of the high-resolution SST field with respect to the coarse SeaFlux/ERA5 SSTs. In order to ensure that the domain-averaged SST correction is zero and that the area-weighted means of the variables in the EUREC⁴A-OA region are conserved, ΔSST takes the following expression

Angle brackets denote the spatial average over the EUREC⁴A-OA region (Ω), namely

Hence, to refine an existing LHF data set we first remap the original SeaFlux/ERA5 U_{10 m} , q_{2 m} , T_{2 m} and SST to the MUR-JPL grid using the closest neighbor technique with Climate Data Operator (CDO) as described in Schulzweida et al. (2019). This is to avoid interpolated values which could potentially affect the results. Then, from these remapped data we compute LHF by means of COARE3.5 (LHF _LR ). Finally, we apply eq. 4 to the remapped SeaFlux/ERA5 variables and introduce them as an input in the COARE3.5 function. This operation provides the refined LHF data set (LHF _HR ).

To validate the downscaling procedure explained above and infer confidence in our results, numerical simulations using the WRF model have been used. Each simulation of a ten member ensemble is run from the 1st of December 2019 to the end of February 2020, forced at the lateral boundaries by ERA5 hourly data. At the lower boundary, the MUR-JPL SST product interpolated on the WRF numerical grid is used (in the following referred to as SST_WRF ). The month of December 2019 is used as a spin-up and is not considered in the following analyses. The validation consists in creating lower resolution data starting from WRF outputs on the original model grid (at 3 km of resolution), that are then refined through the downscaling procedure described above. The obtained LHF _HR field can be then compared with the original LHF _WRF computed using the COARE3.5 formulation directly from WRF variables at 3 km grid spacing. To do so, we have produced lower resolution data, named SeaFlux-like data, by filtering WRF fields U_{10 m} , q_{2 m} , T_{2 m} and SST with a Gaussian filter considering σ = 24km. Such filtered fields are found to have spectral properties that are similar to the original SeaFlux data (not shown). These SeaFlux-like variables correspond to ψ variables in eq. 2, that we thus handle as if they were the original SeaFlux variables: we first compute the associated coupling coefficient and then we apply eq. 4. In the latter step, ΔSST is defined considering SST_HR ≡ SST_WRF .

Figure 1 shows the SeaFlux 2008-2018 DJF climatology in the EUREC⁴A-OA region. Like in most tropical regions, the atmospheric circulation is governed by the trade winds blowing from the open ocean towards South America. Their intensity decreases as they approach the coast ( Figure 1B ). SST and q_{2 m} present the same southwestern-northeastern gradient, both showing larger values towards the coast ( Figure 1A ). An exception to the SST general pattern is the presence of a cold patch close to the South American coast ( Figure 1A ) which likely corresponds to a coastal upwelling system (Acquistapace et al., 2022).

The SST increase towards the coast is a consequence of a western boundary current (the North Brazil Current, NBC) and the anticyclonic eddies (the NBC rings) that this current regularly sheds, which advect northward tropical warm water towards the Caribbean Sea (Johns et al., 1990; Richardson et al., 1994). The saturation water vapor specific humidity is smaller in the open ocean, in line with what expected by the Clausius-Clapeyron relation for evaporation over colder waters. The q_{2 m} and SST fields, together with the decreased wind speeds near the coast ( Figure 1B ), imply lower values of LHF in shelf waters than in the open ocean ( Figure 1B ). Based on these spatial differences, the EUREC⁴A-OA domain is traditionally divided into two sub-regions (Stevens et al., 2021): the NBC eddy corridor region (ECR) and the open ocean region (OOR). They are displayed as the red and grey areas in Figure 1 respectively.

To deepen our understanding in the dynamic and thermodynamic differences between these two regions, Figure 2 represents the histograms of the oceanic and near-surface atmospheric variables analyzed in this paper for the whole EUREC⁴A-OA, for the ECR and for the OOR. Figures 2A, B represent the SST distribution and show that the ECR is characterized on average by warmer SSTs than the OOR.

Figures 2C, D represent the specific humidity deficit (Δq) distribution, defined as the difference between q_* and q_{2 m} . They show that Δq is smaller in the ECR, especially in Figure 2C corresponding to SeaFlux. In turn, the whole EUREC⁴A-OA region and the OOR present similar specific humidity deficit distributions. Accordingly, the U_{10 m} distribution represented in Figures 2E, F shows lower values in the ECR and a wider probability density function (PDF) for the OOR with values ranging from 2 m s^-1 to 12 m s^-1 in both data sets.

In agreement with eq. 1, decreased Δqs and U_10ms lead to lower LHFs in the ECR region. This is shown in Figures 2G, H . Finally, we also check that different LHF estimates provided by three possible ways to compute LHF (eq. 1, COARE3.5 and the LHF provided directly by SeaFlux and ERA5) do not present substantial differences for both, SeaFlux and ERA5 ( Figures 2I, J ). No difference is observed in SeaFlux between the three estimates ( Figure 2I ) since this satellite-based product uses COARE3.0 to derive LHF (Clayson et al., 2014). This algorithm does not differ from COARE3.5 in the turbulent-flux calculation part. For ERA5, we observe slight differences between COARE3.5 and the LHF provided by ERA5 and eq. 1 ( Figure 2J ). These differences are expected as ERA5 does not use COARE3.5 to compute THFs. However, the deviations are less than 10 W m⁻², much smaller than the characteristic LHF differences discussed in this manuscript. Therefore, all the fluxes of this article are hereon calculated with COARE3.5 only. The bulk formula in eq. 1 is also used for theoretical considerations.

This study assesses how fine-scale SST structures influence LHF in the north-west tropical Atlantic Ocean in winter. In order to do that, it is first necessary to understand how the near-surface atmospheric variables are related to SST. We find that this relationship depends on the spatial scale of interest. The analyses of the SeaFlux satellite-based product show that for scales of the order of or smaller than 150 km, q_{2 m} tends to decrease when SST increases whereas U_{10 m} increases by virtue of the DMM mechanism. The change rates are -0.05 g kg^-1 K^-1 and 0.44 m s^-1 K^-1 respectively as shown in Table 2 . However, these results are highly dependent on the data set analyzed. In ERA5 ( Table 2 ), the q_{2 m} -SST coupling coefficient has the opposite sign to that of SeaFlux (0.48 g kg^-1 K^-1) and α_U is 0.44 m s^-1 K^-1.

Our work also highlights the impact, in these coupling coefficients, of the environmental conditions (such as the strength of the prevailing winds) and dynamics which affect their magnitude. The area we studied is characterized by two dynamically different domains: a relatively quiet and homogeneous open ocean and the North Brazil Current ring (eddy) corridor encompassing the South American continental shelf and slope. In the latter, the small-scale U_{10 m} variations with SST are stronger than in the open-ocean domain (0.44 m s^-1 K^-1 versus 0.15 m s^-1 K^-1 as shown in Tables 3 , 4 ). Conversely, the changes in q_{2 m} are higher in the open ocean than in the eddy corridor, possibly as a consequence of the enhanced fine-scale activity in the latter masking the Clausius-Clayperon relationship (-0.07 g kg^-1 K^-1 in the ECR versus 0.34 g kg^-1 K^-1 in the OOR as depicted in Tables 3 , 4 ). This suggests that the environmental conditions such as the presence of an inversion layer at the top of the MABL and the speed of the prevailing winds as well as a different ocean dynamical environment likely affect the sensitivity of the boundary layer properties to the SST. For those reasons, the coupling coefficients and the amplitude of the LHF anomaly obtained from the downscaling procedure are expected to have different magnitudes in different geographical areas and seasons. In particular, it would be interesting to apply the methodology of this study to regions with stronger SST gradients and important mesoscale activity such as the Gulf Stream, the Kuroshio Extension or the Southern Ocean. Apart from ameliorating our understanding of the SST-MABL processes in these regions it could potentially provide a powerful tool to enhance our knowledge on low cloud formation and on air-sea fluxes, with numerous applications in weather forecasting and climate projections.

The interplay described here is therefore only part of a complex series of interactions involving the coupled ocean–atmosphere system which includes two-way interactions between the MABL and the upper ocean. The SST-induced wind disturbances cause perturbations in surface heat fluxes and upper-ocean mixing that are likely to erode SST anomalies that will feed back into the original wind stress perturbations. Moreover, the upwelling associated with SST-induced wind stress curl perturbations will feed back on the ocean, likely altering the SST Small et al. (2008). These two-way feedbacks are intriguing aspects of the coupled system that can significantly enhance our understanding of ocean–atmosphere interactions.

The fine-scale coupling coefficients we derived are at the core of the downscaling algorithm we developed in this work to assess the small-scale SST influence on LHF. When the influence of the ocean small scales is taken into account, the LHF magnitude increases on average of 33% per °C. Further analyses allowed us to distinguish two different mechanisms at play: the thermodynamic and dynamic contributions. The former includes only the increase of LHF with SST associated to the Clausius-Clapeyron dependence of saturation water vapor pressure on SST. The latter is related to the modification of the vertical stratification of the MABL as a consequence of SST anomalies. From SeaFlux we estimate that they represent a LHF increase of 5% and 28% per °C respectively, in agreement with what the LHF bulk formula and previous research predict. Using Surface Wave Instrument Float with Tracking (SWIFT) and Wave Glider observations measured during the ATOMIC campaign as well as satellite data, Iyer et al. (2014) are able to detect differences up to of 30% in LHF across SST gradients (average SST difference of 0.7°C over a distance of 44 km). Our results show that the LHF amplitude is also affected by the geographical heterogeneity of the coupling coefficients. As a result, LHF in the eddy corridor shows the largest variations (around 43% per °C of SST). To provide substantiation to our results, we tested the downscaling approach against an ensemble of one-month long high-resolution regional atmospheric simulations forced at the lower boundary by a high-resolution SST product and at the open boundaries by ERA5. This test shows that the downscaling reduces by a factor of 2 the difference between the simulated and reconstructed LHF fields proving the soundness of the approach.

Our study, however, does not investigate all the processes involved in the near-surface atmospheric response to fine-scale SST structures. For instance, the downscaling technique assumes that all near-surface variables depend exclusively on SST and disregards the nonlinear cross-terms like the increase of air water vapor content due to wind advection. Furthermore, this work does not account for the surface current effect on LHFs estimate. This fact might result in possible inaccuracies, especially in regions characterized by intense mesoscale activity such as the North Brazil current eddy corridor, where surface currents velocity are high (they can exceed 1 m s^-1) (Seo, 2017; Renault et al., 2019; Zhang et al., 2021). To test how this effect weights in in LHF, we use data from the two NOAA-funded Saildrones of the EUREC⁴A-OA/ATOMIC experiment (Stevens et al., 2021): SD1063 and SD1064. These wind-propelled uncrewed surface vehicles were equipped with different sensors collecting measurements at the air-sea interface (e.g., SST, sea-surface salinity, winds, air temperature…) at 1-minute temporal resolution thereby offering an unprecedented view of the upper ocean and the air-sea interface. They were also equipped with an Acoustic Doppler Current Profiler (ADCP) on the keel, allowing to monitor the upper ocean currents. SD1063 sampled parts of the eddy corridor while SD1064 navigated through the open-ocean sector of the EUREC⁴A-OA domain. Using these in-situ measurements, we find that the intense surface currents at the edge of mesoscale eddies induce a change of almost 5% to 15% in LHF (not shown). This suggests that other data sources are needed in order to further explore the effects of fine-scale SST structures on air-sea fluxes. The high temporal resolution of in-situ data from the EUREC⁴A-OA/ATOMIC campaigns (Quinn et al., 2021; Stevens et al., 2021) are likely to provide better estimates of the impacts of SST gradients on air-sea fluxes that go beyond the coarse (with respect to the ocean small-scale) spatial grid of available satellite observations.

The data sets analyzed for this study can be found in the following open-access repositories.

All the authors contributed to conception and design of the study. SS and PF organized the SeaFlux, ERA5 and MUR-JPL databases and performed the statistical analysis with them. In turn, AM, MB, FD, and CP organized the WRF model database and performed the statistical analysis with it. PF wrote the first draft of the manuscript and the sections which did not involve analysis with the WRF model. AM, MB, FD, and CP wrote sections of the manuscript related to the WRF data description and the downscaling algorithm validation. The separation technique between the thermodynamic and dynamic contributions of LHF sensitivity to SST was performed by AM, MB, and CP. All authors contributed to the article and approved the submitted version.