On the importance of the atmospheric coupling to the small-scale ocean in the modulation of latent heat flux

In this study, ocean and atmosphere satellite observations, an atmospheric reanalysis and a set of regional numerical simulations of the lower atmosphere are used to assess the coupling between the sea-surface temperature (SST) and the marine atmospheric boundary layer (MABL) as well as the latent heat flux (LHF) sensitivity to SST in the north-west tropical Atlantic Ocean. The results suggest that the SST-MABL coupling depends on the spatial scale of interest. At scales larger than the ocean mesoscale (larger than 150 km), negative correlations are observed between near-surface wind speed (U1 0m ) and SST and positive correlations between near-surface specific humidity (q 2m ) and SST. However, when smaller scales (1 – 150 km, i.e., encompassing the ocean mesoscale and a portion of the submesoscale) are considered, U10 m -SST correlate inversely and the q 2m -SST relation significantly differs from what is expected using the Clausius-Clapeyron equation. This is interpreted in terms of an active ocean modifying the near-surface atmospheric state, driving convection, mixing and entrainment of air from the free troposphere into the MABL. The estimated values of the ocean-atmosphere coupling at the ocean small-scale are then used to develop a linear and SST-based downscaling method aiming to include and further investigate the impact of these fine-scale SST features into an available low-resolution latent heat flux (LHF) data set. The results show that they induce a significant increase of LHF (30% to 40% per °C of SST). We identify two mechanisms causing such a large increase of LHF: (1) the thermodynamic contribution that only includes the increase in LHF with larger SSTs associated with the Clausius-Clapeyron dependence of saturating water vapor pressure on SST and (2) the dynamical contribution related to the change in vertical stratification of the MABL as a consequence of SST anomalies. Using different downscaling setups, we conclude that largest contribution comes from the dynamic mode (28% against 5% for the thermodynamic mode). To validate our approach and results, we have implemented a set of high-resolution WRF numerical simulations forced by high-resolution satellite SST that we have analyzed in terms of LHF using the same algorithm. The LHF estimate biases are reduced by a factor of 2 when the downscaling is applied, providing confidence in our results.


Introduction
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 .
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 Boâs 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 stateof-the-art bulk parametrization outputs (Shao et al., 2019). Cloudresolving 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 4 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 4 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 airsea interface has been achieved. This was made possible by additional initiatives such as the European EUREC 4 A-OA (Karstensen et al., 2020;Speich et al., 2021, EUREC 4 A-OA) and the American ATOMIC (Quinn et al., 2021, Atlantic Tradewind Ocean-atmosphere Mesoscale Interaction Campaign) ones. In particular, EUREC 4 A-OA constitutes the ocean component of EUREC 4 A. It investigates heat, momentum, freshwater and CO 2 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 4 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 4 A-OA, this study focuses on the 5°-17°N, 60°-51°W box, hereby referred to as the EUREC 4 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.

SeaFlux
The SeaFlux project (Clayson et al., 2014) is one of the multiple ongoing efforts to develop satellite-based estimates of the ocean surface THFs and associated near-surface properties. These efforts principally rely on microwave imager observations from which nearsurface wind speed, specific humidity and temperature are estimated (Bourassa et al., 2010). The SeaFlux version used in this paper is called SeaFluxV3, which makes use of a nonlinear neural network technique to estimate near surface air properties from microwave radiances . Moreover, diurnally varying SSTs are obtained superimposing a sinusoidal diurnal cycle onto the foundation sea surface temperature provided by the NOAA Optimally Interpolated SST (Reynolds et al., 2007). The reference temperature is taken as the pre-dawn SST and the lag between noon and the maximum warming is estimated to be 0.28 times the length of the day (Clayson et al., 2014). Finally, to compute surface turbulent fluxes from near-surface variables, SeaFlux makes use of a neural network version of COARE3.0 algorithm (Fairall et al., 2003) which has been optimized for speed of use of processing. No additional parametrizations such as those taking into account the effects of waves or sea spray are considered.
This data set has a spatial resolution of 25 km × 25 km and a time resolution of 1 hour. The variables we use in the present study are 2 m specific humidity (q 2m ), 10 m wind speed (U 10m ), 2 m air temperature (T 2m ), SST and LHF. Only daily averages of all the variables are considered in this study.

ERA5
The ERA5 reanalysis (Hersbach et al., 2020) embodies a detailed record of the global atmosphere, land surface and ocean waves from 1950 onward. It does not include the coupling with an ocean model, but it is forced at the lower boundary by a prescribed SST and seaice concentration. In particular, it exploits various flavors of the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) data set as well as the Climate Change Initiative (ESA-CCI) SST v1.1 (Merchant et al., 2014) up to 2007. The Operational Sea Surface Temperature and Ice Analysis (OSTIA) is considered for the modern period (Donlon et al., 2012).
ERA5 has a spatial resolution of 0.25°× 0.25°and a time resolution of 1 hour. The variables used are 2 m dew point temperature, SST, T 2m , 10 m wind components, sea-level pressure (SLP) and LHF. U 10m is computed taking the norm of the horizontal wind components and q 2m is obtained from 2m dew point, T 2m and SLP as in Buck (1981). Only daily averages of all the variables are considered in this study.

MUR-JPL
This data set provides high-resolution SST values distributed over a global 0.01°× 0.01°grid. SST values from version 4 Multiscale Ultrahigh Resolution (MUR) Level 4 analysis (Chin et al., 2017) are based upon nighttime observations from several instruments including the NASA Advanced Microwave Scanning Radiometer-EOS (AMSR-E), the JAXA Advanced Microwave Scanning Radiometer 2 on GCOM-W1, the Moderate Resolution Imaging Spectroradiometers (MODIS) on the NASA Aqua and Terra platforms, the US Navy microwave WindSat radiometer, the Advanced Very High Resolution Radiometer (AVHRR) on several NOAA satellites, and in-situ SST observations from the NOAA iQuam project. Daily outputs are obtained with a multi-scale variational approach that combines all the available observations.

WRF
The Weather Research and Forecasting (WRF) model V4.1.5 (Skamarock et al., 2019) is used with its Advanced Research WRF (ARW) core. It solves the non-hydrostatic fully compressible Euler equations on an Arakawa-C grid with hybrid vertical coordinates. Two double-way nested domains with a grid spacing of 9 and 3 km and 75 vertical levels (25 of them are in the lowest 2 km) are used in our study. The innermost domain covers the area used for the SeaFlux and ERA5 analyses, namely the EUREC 4 A-OA region. Conversely, the outermost domain spans from 4°S to 21°N and from 67°W to 31°W. The simulation used here is a part of a set of experiments designed specifically for the EUREC 4 A-OA region and has been submitted with some details about the configuration. Briefly, boundary conditions, updated every hour, are provided by ERA5 for the outer domain. The cloud microphysics parametrization is given by the Thompson-Eidhammer scheme (Thompson and Eidhammer, 2014), the planetary boundary layer (PBL) is parameterized using the Yonsei University PBL representation (Hong et al., 2004), the radiation parametrization follows the Rapid Radiative Transfer Model (Mlawer et al., 1997;Iacono et al., 2008) and cumuli are represented with the Kain-Fritsch algorithm (Kain, 2004) for the outer domain which includes subgrid-scale cloud feedback to radiation. Finally, the SST is prescribed daily by means of the MUR-JPL SST (Chin et al., 2017).

Methodology
THFs are often computed using bulk formulae. In particular, a widely-used expression for LHF (Yu, 2009) where r 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 10m , 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 10m , q 2m , T 2m and SST are filtered with an isotropic Gaussian filter as described in the Appendix. Hence, a given variable y is decomposed as with y being the filtered large-scale component and y ′ the residual field, which contains small-scale information. The standard deviation of the Gaussian filter s represents the threshold imposed to separate the small and the large-scale during the filtering (see the Appendix). Unless otherwise stated s = 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 airocean coupling. We indeed compute, through a linear regression, the coupling coefficients between SST' and other variables such as U' 10m , q' 2m and T' 2m . The resulting coupling coefficients are denoted TABLE 1 COARE3.5 input parameters for SeaFlux and ERA5: Sea-level pressure (SLP), downward shortwave radiation (SW), downward longwave radiation (LW), wind sensor height (z u ), bulk temperature sensor height (z t ), humidity sensor height (z q ), mean latitude of the EUREC 4 A-OA region and MABL height.

Variable
Value Mean latitude (°N) 12.5 MABL height (m) 600 as a U , a q and a 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 4 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 (y) we build a new data set of the same variable including fine-scale SST features (y HR ) as where DSST 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 4 A-OA region are conserved, DSST takes the following expression Angle brackets denote the spatial average over the EUREC 4 A-OA region (W), namely Hence, to refine an existing LHF data set we first remap the original SeaFlux/ERA5 U 10m , q 2m , T 2m 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 10m , q 2m , T 2m and SST with a Gaussian filter considering s = 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 y 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, DSST is defined considering SST HR ≡ SST WRF . Figure 1 shows the SeaFlux 2008-2018 DJF climatology in the EUREC 4 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 2m 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).

Results
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 2m 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 4 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 4 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 (Dq) distribution, defined as the difference between q * and q 2m . They show that Dq is smaller in the ECR, especially in Figure 2C corresponding to SeaFlux. In turn, the whole EUREC 4 A-OA region and the OOR present similar specific humidity deficit distributions. Accordingly, the U 10m 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 Dqs and U 10m s 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 −2 , 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.

SST-MABL coupling at different spatial scales
Given the spatial heterogeneity of the EUREC 4 A-OA region shown in Figures 1, 2, point-wise regressions are performed between SST and q 2m and U 10m . They are shown in Figure 3. White contours in Figures 3A, B show the regression coefficient between the smoothed SST and q 2m fields for SeaFlux and ERA5 respectively. In both cases, the large-scale near surface specific humidity increases with SST in all the locations of the EUREC 4 A-OA region in agreement with the fact that warmer air needs more water vapor to reach saturation. The estimated coupling coefficients range between 0.84 g kg -1 K -1 and 1.27 g kg -1 K -1 , slightly less than the Clausius-Clapeyron scaling a q*~1 .3 g kg -1 K -1 . Such a scaling is obtained with a linearization of the Clausius-Clapeyron expression of Buck (1981) (as used in COARE3.5) around the reference SST value of 300 K, representative of the region (see the Appendix for the full derivation). Physically, this corresponds to the fact that at large scales there is enough evaporation so that the near-surface specific humidity is close to saturation.
Conversely, the coupling coefficients of the SeaFlux SST-q 2m residual fields indicate that, at the small-scale, the Clausius-Clapeyron scaling is not respected (Figures 3A, B). The SeaFlux residual coupling coefficients range from -0.4 to -0.1 g kg -1 K -1 and the ERA5 ones lie between 0.1 and 0.4 g kg -1 K -1 , in both cases significantly smaller than the 1.3 g kg -1 K -1 reference value. This fact is interpreted in terms of a more complex atmospheric response to the small-scale oceanic forcing. In fact, warmer SSTs are associated with a more unstable MABL. This configuration enhances vertical mixing within the MABL, entrainment at the MABL top and a subsequent MABL thickening. Hence, at least two mechanisms contribute to the decrease of a q : (1) the thickening of the MABL that dilutes the air water content, and (2) the mixing at the top of the MABL which entrains dryer air from the free atmosphere (Neggers et al., 2006;Acquistapace et al., 2022;Albright et al., 2022).
Note that the SeaFlux q 2m coupling coefficients are negative over most of the EUREC 4 A-OA domain ( Figure 3A) whereas the ERA5 ones remain slightly positive ( Figure 3B). This fact might suggest that the MABL dynamical mechanisms which reduce the surface specific humidity are weaker in ERA5 than in SeaFlux. It is likely that the numerical parametrizations used in the ERA5 model are not capable to correctly reproduce the MABL dynamics. In particular, the schemes that parameterize the surface fluxes and the MABL turbulence might play a major role (Perlin et al., 2014). Elucidating the specific causes of this difference requires further research which is not presented here as it goes beyond the scope of this study. What matters for our purposes is that the coupling EUREC 4 A-OA region (5°-17°N, 60°-51°W) DJF climatology for SeaFlux. The shading in (A) represents the climatological SST and the contours stand for near-surface specific humidity (in g kg -1 ). In (B) the shading represents air-sea latent heat flux and the contours illustrate the mean values of near-surface wind speed (in m s -1 , white color is used close to the coast to facilitate visualization). In both maps the gray box delimits the OOR (14°-16°N, 58°-52°W) and the red box represents the ECR. Daily values between 2008 to 2018 are considered to compute the climatologies.
coefficients of the specific humidity from both SeaFlux and ERA5 residuals are different from the Clausius-Clapeyron scaling, consistent with the fact that air moving over small SST structures does not have the time to reach an equilibrium condition with the sea underneath. The large-scale U 10m and SST linear regression slopes are shown in black contours in Figures 3C, D for SeaFlux and ERA5 respectively. They are negative over almost the whole EUREC 4 A-OA region. They reach their minimum over the NBC retroflection, perhaps as a result of the increase in SST and the weakening of surface winds as we approach the equator. The negative coupling coefficients are often interpreted as a consequence of the ocean responding passively to the atmospheric forcing. Higher nearsurface wind speeds are responsible for increased LHF and enhanced mixing in the upper ocean. These two factors result in a decrease of SST.
Instead, positive coupling coefficients for the residual U 10m and SST fields dominate over the EUREC 4 A-OA region, especially in the ECR (shading in Figures 3C, D). In the attempt to explain how small-scale SST features affect surface winds, two mechanisms have been previously identified: the pressure adjustment (PA) mechanism (Lindzen and Nigam, 1987)  Histograms of the different variables analyzed in this paper for SeaFlux (left column) and ERA5 (right column). (A, B) represent the SST; (C, D) depict the 2 m specific humidity deficit defined as the difference between the saturation specific humidity and q 2m ; (E, F) illustrate the U 10m distribution and (G, H) display LHF computed using COARE3.5. In all the plots hitherto mentioned the green curve is obtained considering the whole EUREC 4 A-OA box (5°-17°N, 60°-51°W), and the blue and red histograms are calculated over the OOR and ECR regions respectively as defined in the main text and as shown in Figure 1. (I, J) represent the LHF given directly in the data sets (green), computed using eq. 1 (blue) and using COARE3.5 (red). All the histograms are constructed using daily averages of SeaFlux and ERA5 data in the 2008-2018 DJF season.
momentum mixing (DMM) Wallace et al., 1989). In the PA mechanism the thermal expansion of air over warm SST patches generates pressure anomalies that drive a secondary circulation characterized by surface wind convergence at the center of the warm anomaly. Similarly, surface wind divergence occurs over cold anomalies. This implies weak surface wind velocities both at SST maxima and minima, and therefore no correlation SST-U 10m . In turn, the DMM mechanism expresses that, when an air parcel crosses a SST front going from cold to warm water, its buoyancy increases and reduces the static stability of the air column. This enhances a downward flux of horizontal momentum from the top of the MABL to the surface, triggering the entrainment of higher values of wind speed into the MABL thereby increasing surface wind speed on the warm side of the front, this results in a positive SST-U 10m correlation. The positive and statistically significant a U s shown in Figures 3C, D, especially in the ECR, suggest the DMM mechanism is here at play and it dominates the small-scale SST-U 10m interactions. A more detailed discussion on the two mechanisms, involving different metrics computed in the EUREC 4 A-OA region, can be found in the Appendix. Spatial distribution of the and SST-q 2m and SST-U 10m linear regression slopes for SeaFlux (left column panels) and ERA5 (right column panels). (A, B) represent the SST-q 2m linear regression slopes for the smoothed (contours, in g kg −1 K −1 and in white to facilitate visualization) and residual fields (shading). (C, D) represent the SST-U 10m linear regression slopes for the smoothed (contours, in m s −1 K −1 ) and residual fields (shading). In all the cases, the hatching marks the slopes from the residual fields which are statistically significant at a 99% confidence interval after a two-sided t-test. All the panels are obtained considering daily data during the 2008-2018 DJF season. The s value introduced in the Gaussian filter equals 150 km. In all the maps the OOR and the ECR are delimited by the northern and coastal boxes respectively. So far, the spatial distribution of the coupling coefficients has been assessed for a fixed value of s = 150 km. To obtain a deeper insight into the dependence of the SST-MABL coupling on the spatial scale, we consider a wider range of s values. In particular, we compute a U and a q for s ranging from 100 km to 1000 km ( Figure 4). The coefficients computed for ss smaller than 100 km (the effective resolution of the data sets according to the autocorrelation analysis of the Appendix) are not presented in Figure 4 as they do not contain anything but random noise. The sub-sampling spacing to compute the coupling coefficients presented in Figure 4 is taken as 150 km regardless of the s value. The reasons for this choice are explained in the Appendix. Figure 4A shows that a q monotonically grows with increasing s for the two data sets. In ERA5, it tends to the Clausius-Clapeyron scaling value of 1.3 g kg -1 K -1 . Apart from the difference in value (the SeaFlux coupling coefficient is always 0.5 g kg -1 K -1 -0.7 g kg -1 K -1 smaller than ERA5), all the data sets provide low values of the coupling coefficients for small values of s. This suggests that, at small-scales, specific humidity does not adjust to SST anomalies with values close to saturation as relative humidity changes and advection probably mask the direct relationship expected from Clausius-Clapeyron. Instead, the dynamical MABL response to the small-scale SST structures reduces the specific humidity dependence on SST, because of the modulation of the MABL height and its mixing. In turn, the SST-U 10m coupling coefficients for SeaFlux and ERA5 are generally positive, as expected from the DMM mechanism ( Figure 4B). In addition, as s increases, the coupling coefficients of SeaFlux and ERA5 decrease as larger-scale features are retained in the residual fields, in agreement with the change of atmosphere-ocean coupling between large and small scales discussed in section 1.
Note that both SeaFlux SST-q 2m and SST-U 10m coupling coefficients present a tendency change at around s = 200 km. As we approach s = 200 km from larger s values we observe a change in sign from positive to negative in a q ( Figure 4A) and the maximum value of the a U -s scatter plot ( Figure 4B). Such behaviors are consistent with all the discussion of the preceding paragraphs regarding the small-scale SST-near-surface atmosphere linkages as 200 km is roughly the limit between the ocean mesoscale and the large-scale. This fact also supports our choice of 150 km as the separation scale to define the residual fields.
The estimated coupling coefficients will be used in the following to downscale the impact of the small-scale SST on the LHF. Even though they have been proven to be sensitive to the spatial scale considered and the specific location within the EUREC 4 A-OA region, for simplicity, only the coupling coefficients in Figure 4 for s = 150 km have been used in the rest of the work and are shown in Table 2. However, for the sake of completeness, the coupling coefficients in the ECR and OOR as defined in Figure 1 for s = 150 km are also provided in Tables 3, 4 respectively. All the calculations in next section are also performed with these coupling coefficients as well in order to give an estimate of the variability of the results as a function of the spread in the coupling coefficients.

LHF modulation by air-sea coupling and sensitivity to SST
Figures 5A, B represent the averaged 2008-2018 DJF difference between LHF HR and LHF LR for SeaFlux and ERA5 respectively. Recall that for SeaFlux and ERA5, LHF HR represents the COARE3.5 latent heat flux output computed from the U 10m , q 2m , T 2m and SST given by eq. 4 and that LHF LR is the COARE3.5 latent heat flux output obtained from the remapped U 10m , q 2m , T 2m and SST. These two LHFs show differences up to 10 W m -2 in coastal areas of the ECR. Weaker LHF values prevail by moving towards the OOR. The quantities in brackets over Figures 5A, B represent the areaweighted means of the LHF differences. They are close to zero for both datasets. This means that, away from the locations where DSST is large (i.e. the ECR), the O (DSST 2 ) effects in the LHF controlling variables (U 10m , q 2m , SST and T 2m ) do not produce large LHF variations so that LHF evolves almost linearly with DSST. Accordingly, the spatial pattern of the LHF differences coincides with the size of the DSSTs shown in Figures 5C, D for SeaFlux and ERA5 respectively.  The differences of 10 W m -2 shown in Figures 5A, B only represent around 5 -7% of the mean LHF in the EUREC 4 A-OA region (shading in Figure 1B). These deviations are reduced when averaged over the 10 DJF seasons because the characteristic lifetime of the fine-scale ocean features whose effects on LHF we aim to analyze ranges from days to months. As a consequence, it is interesting to look at daily snapshots of the LHF bias and DSST in order to study the response of LHF to these fine-scale SST structures. Let us then consider Figure 6 which depicts the same fields as Figure 5 but for the 1 st of December 2008. Biases up to 60 W m -2 in absolute value are observed, (Figures 6A, B) especially in the ECR which are equivalent to 30% to 40% of the climatological LHF (shading Figure 1B), responding to DSSTs not exceeding 1°C in absolute value ( Figures 6C, D). Again, the spatial pattern of the LHF differences resembles to that of the SST correction suggesting a linear relationship between the two quantities. We come back to this point in the next paragraphs.
So far, the downscaling has been applied using the coupling coefficients computed accounting for the whole EUREC 4 A-OA region. However, the coupling coefficients present strong spatial variations within the EUREC 4 A-OA region (Figure 3) which have an impact in the resulting downscaled fluxes. To provide an estimate on how large this effect is, we downscale the LHF using the coefficients in Tables 3, 4 corresponding to the ECR and OOR respectively. Strong deviations of the order of 40% to 50% are observed in both cases for particular days of the 2008-2018 DJF period, evincing the importance of maintaining the distinction between sub-regions in the downscaling.
Given the large LHF variations with SST, we try here to separate the different mechanisms by which the fine-scale SST structures may influence LHF. This concept is hereafter referred to as LHF sensitivity to SST. We group these mechanisms into two main categories: the thermodynamic contribution and the dynamic contribution. The thermodynamic contribution includes only the increased LHF at larger SST due to the Clausius-Clapeyron dependence of saturation water vapor pressure to SST, while the surface relative humidity and the air-sea temperature difference are considered to be fixed. The dynamic contribution relates to the MABL response to increased SST, which reduces air-column stability, increases the vertical mixing in the lower troposphere, and favors entrainment at the top of the MABL, resulting in increased surface winds and specific humidity deficit, both of which enhance the LHF from the ocean into the atmosphere.
In order to theoretically quantify the importance of each contribution, let us consider the bulk formula for LHF (eq. 1). We choose the most frequent values of the corresponding histograms in Figure 2 as representative quantities for the EUREC 4 A-OA region: U 10m = 9m s -1 , q * = 21.7g kg -1 (saturation specific humidity at 27°C) and q 2m = 15.5g kg -1 , C e = 1.6·10 -3 . This leads to LHF ≃ 160 W m -2 using eq. 1. Consider now that SST is increased by 1°C and compute first the thermodynamic change of LHF: the only term that changes is (q * − q 2m ) which, assuming a constant relative humidity, increases by ∂q * /∂SST(1−q 2m /q * )≃0.31 g kg -1 K -1 . Considering that (q * − q 2m ) is about 5.5 g kg -1 , this results in an increase in LHF of about 5% (=0.31/5.5), per the 1°C increase in SST.
Using the SeaFlux coupling coefficients shown in Table 2, we compute the effect of the dynamic response on LHF: U 10m increases by 0.44 m s -1 K -1 , which means a relative increase of 4.9% (=0.44/9), (q * − q 2m ) increases by about 1.3 g kg -1 (we consider here q 2m vs SST slope equal to 0, 1.3 g kg -1 is the increase in q * when SST changes from 26°C to 27°C) -relative increase of (q * − q 2m ) is about 24% per°C (=1.3/5.5). In total, the sensitivity of LHF to SST through the dynamic response is about 28% per°C, a factor of about 4-5 larger than the thermodynamic response. Putting all together, the sensitivity of LHF to SST is about 33% per°C.
We perform the same calculation with the SeaFlux coupling coefficients shown in Tables 3, 4 in order to have an estimate of the spatial variability of the LHF sensitivity to SST within the EUREC 4 A-OA region. According to Figure 2, the ECR has a typical SST of 27.1°C, Dq is about 4.1 g kg -1 with q *~2 2 g kg -1 and q 2m~1 8 g kg -1 . Moreover, U 10m = 8m s -1 and a U~0 .44 m s -1 K -1 . These values produce a thermodynamic and dynamic contributions of around 6% and 37% of LHF change per°C of SST respectively. These two contributions together represent a change of 43% of LHF change per°C. Conversely, for the OOR we have SST~of 26.8°C, Dq is about 5.8 g kg -1 with q *~2 1.5 g kg -1 and q 2m~1 6 g kg -1 . Moreover, U 10m = 9m s -1 and a U~0 .15 m s -1 K -1 . These values result in a thermodynamic contribution of around 6% and a dynamic contribution of 24% producing a total sensitivity of LHF to SST of the order of 30%. Therefore, LHF is not impacted by SST variations in the same way in the entire EUREC 4 A-OA region. The difference comes mainly from the fact that the MABL response to fine-scale SST structures is weaker in the OOR than in the ECR.

SeaFlux ERA5
a T (unitless) 0.60 a T (unitless) 0.53 We also test the results from the last two paragraphs using SeaFlux and ERA5 data. To compute the total LHF sensitivity to SST we perform a linear regression between the relative change in LHF (100· (LHF HR − LHF LR )/LHF LR ) and DSST. The results for SeaFlux and ERA5 are shown in Figure 7A in blue and red respectively. The estimated LHF change per°C for SeaFlux (33.8% per°C) agrees with the theoretical value mentioned above whereas the ERA5 LHF sensitivity to SST is weaker (26.6% per°C). The shading in Figure 7A indicates the range of different LHF variations obtained between the OOR (lower limit) and the ECR (upper limit) for SeaFlux (ERA5 results are not included as they do not modify the analysis and conclusions). It aims to provide an estimate of the uncertainty in the linear regressions associated to the geographical heterogeneity of the MABL response to fine-scale SST structures. For small values of DSST, the LHF variations are similar within the EUREC 4 A-OA region. However, for |DSST| larger than 0.5°C we obtain LHF variations ranging from ±20% -±40% in the OOR and ECR respectively. Figure 7A also allows us to identify that for SST corrections with absolute values lower than 0.5°C, the LHF response is almost linear despite COARE3.5 being a non-linear algorithm. Beyond the ± 0.5°C threshold, LHF changes become stable around ±20% -±30%.
In order to isolate the thermodynamic contribution in the LHF sensitivity, a new LHF data set (LHF therm ) is produced with COARE3.5. In this case, only SST and T 2m are modified by adding the SST correction (the coupling coefficients are not involved here). This is to keep the vertical MABL temperature gradient, q 2m and U 10m as constant values. Again, we perform a linear regression between the relative change in LHF and DSST and study the value of the slopes ( Figure 7B). In agreement with the theoretical result, the thermodynamic contribution accounts for the 5.2% per°C of SST of LHF change in SeaFlux. Accordingly, ERA5 provides a slightly weaker sensitivity of 4.7% per°C of SST. In Figure 7B we also appreciate the boundaries of the linear regime and the importance of treating separately the different sub-regions of the EUREC 4 A-OA domain to accurately determine the LHF variations associated to finescale SST structures, especially for large values of |DSST|.

Validation
We validate our approach by means of the original highresolution LHF obtained from the original simulated WRF variables using the COARE3.5 algorithm ( Figure 8A). This LHF dataset, which we denote as LHF WRF , is physically consistent with the fully compressible Euler equations and, thus, can be considered a reliable source of high-resolution information.
To test the robustness of the downscaling approach, we look at the difference between LHF WRF and LHF LR,HR,therm . To estimate the downscaled LHFs we use the coupling coefficients computed from LHF LR that do not differ significantly from Seaflux (not shown). Results in Figure 8 show that LHF HR obtained from the downscaling best represents LHF WRF both in terms of spatial patterns ( Figure 8D) and in terms of spatial-temporal statistics ( Figure 8B). Figure 8B displays the PDFs of differences of LHF LR (blue curve), LHF therm (green curve) and LHF HR (orange curve) with LHF WRF ; the latter is the narrowest PDF among the three showing that the downscaling approach we have implemented to isolate the small-scale ocean influence on the atmosphere reduces (approximately by a factor 2) the biases between LHF datasets.
Concerning the spatial structure of LHFs, Figures 8D-F display the differences LHF HR -LHF WRF , LHF LR -LHF WRF and LHF therm -LHF WRF respectively. They show that LHF HR accurately represents the reference LHF both in magnitude and in the spatial distribution of the bias. The spatial patterns in Figures 8C, E, which are similar to each other, are determined by DSST ( Figure 8C) confirming the linear relationship of the downscaling as in the case of Figure 7A. However, the magnitude of the correction provided for the downscaling applied to the WRF simulated fields is smaller than in the case of SeaFlux. This is due to the fact that the range of DSST correction ( Figure 8C) is much smaller for WRF than for SeaFlux in Figure 7A.
It is also interesting to note that the correction DSST in Figure 8C differs from Figures 5C, D also in terms of spatial patterns. The former presents only small-scales correction, while in the latter largescales structures emerge. This might be due to the contribution to the correction in Figures 5C, D from two different products (SeaFlux and MUR-JPL) while Figure 8C is obtained only with MUR-JPL. The inconsistencies and biases between these two different data sets do not allow us to directly compare the spatial SST correction patterns and might also explain the noisy pattern shown in Figure 5.

Summary and conclusions
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 2m tends to decrease when SST increases whereas U 10m 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 2m -SST coupling coefficient has the opposite sign to that of SeaFlux (0.48 g kg -1 K -1 ) and a 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 10m 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 2m 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 s i gn ifi c an t ly e n h an c e ou r u n d e r s t an d i n g o f oc e anatmosphere 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 Validation with WRF data. In the figure, LHFs are all computed by means of the COARE3.5; LHF WRF is the LHF computed from WRF atmospheric and oceanic outputs, LHF LR represents the LHF computed from upscaled WRF data (using a Gaussian filter with s = 24 km) and LHF HR stands for the LHF data set computed after having applied the downscaling technique to q 2m , U 10m , T 2m and SST. LHF therm represents the only thermodynamic contribution to LHF, obtained from the correction of T 2m and SST with DSST. (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 4 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 4 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 4 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.

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

Author contributions
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.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Appendix : Statistical definitions, methods and further considerations
Data filtering A given variable y can be convoluted over the sea with a Gaussian filter defined as: G(x, y) = exp ( − x 2 +y 2 2s 2 ) for (x 2 + y 2 ) 1 2 < 3s 0 otherwise In this way, the filtered variable is with a normalization factor The standard deviation of the Gaussian filter s is tightly linked to the cutoff length imposed by the filter itself. Using classical Fourier transform arguments (Arfken and Weber, 2005) it can be shown that the cutoff wavenumber (i.e. the maximum resolved length of the filter) is inversely proportional to s. The reader is referred to Meroni et al., (2018) for further details.

Autocorrelation length
The procedure followed to compute the autocorrelation length is detailed in Meroni et al. (2020) and briefly outlined in this section. Given a generic field q(j,q) in spherical coordinates (j is the longitude and q is the latitude), the bidimensional autocorrelation function is defined as Here {x, h} is a set of standard local Cartesian coordinates (x is positive eastward and h is positive northward). |W|, s 2 q and q represent, respectively, the area of the region of interest, the spatial variance and the area-weighted mean of the field, namely: q = 1 W j j ðð W dqdjR 2 cos qq(j, q) : The isotropic autocorrelation function A q (r) is found on a set of local polar coordinates {r, d}. r is the radial distance from the origin and d is the counterclockwise angle from the positive x. A q (r) is obtained by averaging A q (x,h) over full circles, namely A q (r) = 1 2p Z 2p 0 dd A q (r cos d , r sin d ) : Once the isotropic autocorrelation function is known, the autocorrelation length L is obtained by solving A q (L) = e −1 : The autocorrelation length is computed on the daily residual fields. Its temporal average over the 2008-2018 DJF season is taken as an estimate for the sub-sampling distance in the computation of the coupling coefficients. Note that the autocorrelation length is linked to the s value chosen in the filter. Ideally, the autocorrelation length ought to be computed for the residual fields given by each s in order to obtain a statistically consistent subsampling spacing. However, this procedure does not leave enough grid-points to perform the regressions and obtain robust results of the coupling coefficients for s> 600 km as the EUREC 4 A-OA region is roughly 10°wide in latitude and longitude. Therefore, for all the residual fields, the autocorrelation length is chosen to be 150 km. Only the grid-points which are separated by a distance corresponding to the autocorrelation length are chosen to perform the linear regressions. This procedure enables us to ensure that all the values used in the regressions are statistically independent.

Linearization of the Clausius-Clapeyron equation
Consider the integrated Clausius-Clapeyron equation written in terms of the saturation water vapor pressure as used in COARE3.5: e * = e * 0 exp ½ B · T C + T , where e *0 = 6.1121 hPa, B = 17.502°C -1 and C = 240.97°C three constants of integration. T stands for temperature and is expressed in degrees Celsius. For simplicity no salinity and pressure-related corrections are considered in the following steps as the results do not change much (not shown). It is standard practice to consider SST as the controlling temperature for the computation of the surface latent heat flux. Let us write it as the sum of a large-scale field SST and the residuals SST ′ , namely SST = SST + SST 0 . The saturation water vapor pressure is, then, e * = e * 0 exp ½ B · (SST + SST 0 ) C + SST + SST 0 which corresponds to e * (x) = e * 0 exp ½ B(1 + x) being x the ratio between SST ′ and SST. This can be expanded around zero as: With e * (0) = e * 0 exp ½ B · SST C + SST , that is the large scale saturation water vapor pressure field, and de * dx j x=0 = e * (0) BC SST½ C SST+1 : Thus, the saturation water vapor pressure is e * (SST 0 ) = e * (0)½1 + B · C ( C SST + 1) 2 with the first term that depends on SST only, and the second term, the spatial residual, that is linearly dependent on the SST residuals, SST ′ , with a slope controlled by the large scale SST, i.e. SST. Using the approximation that the atmospheric pressure is a smooth field, denoted with p, the saturation specific humidity is readily written as q * (SST 0 ) ≈ ϵ e * (SST 0 ) p − e * = ϵ e * p − e * + ϵ e * p − e * B · C ( C SST + 1) 2 SST 0 SST 2 (23) with ϵ = 0.622 the ratio between the dry air and the water vapor gas constants. Therefore one can express the smoothed saturation specific humidity q * such that: Taking the derivative of the second term of eq. 23 with respect to SST ′ we obtain the coupling coefficient of the saturation water vapor specific humidity: Thus, considering SST ∼ 26:8 ∘ C we get q * ∼ 23 Â 10 −3 kg kg -1 and an easy calculation introducing all the values provided in this section in eq. 25 provides a value of ∂q ′ * /∂SST ′ of 1.3 g kg -1 K -1 . This reference value is widely used throughout the manuscript.

Role of SST structures on MABL dynamics
Binned scatter plots of (A) along-wind SST gradient (∂SST/∂r) and along-wind wind divergence (∂ur/∂r), and (B) across-wind SST Laplacian (∂ 2 SST/∂s 2 ) and across-wind wind divergence (∂us/∂s). Transparent symbols indicate that less than 100 values were considered in the corresponding bin. In the lower part of the figure, the histogram of the forcing SST fields are shown. The legend contains the linear regression slope of the scatter plots denoted as a and its p-value in parentheses.