Air-sea turbulent heat flux feedback over mesoscale eddies

so-called


Introduction
The turbulent heat flux feedback (THFF, in  hereafter) is a critical parameter, which measures the change in the net air-sea turbulent heat flux in response to a 1 K change in sea surface temperature (SST).It is a powerful tool to quantify the rate of dampening of SST anomalies.THFF can vary seasonally (largest in winter), geographically and with ocean spatial scale.Early studies estimate THFF at ∼20 2 1 W m K E   for basin-scale mid-latitude SST anomalies, which, to first order, respond passively to atmospheric forcing (Bretherton, 1982;Frankignoul, 1985;Frankignoul et al., 1998Frankignoul et al., , 2004;;Small et al., 2020).More recent studies estimate that THFF increases to 40 W m K E   in the Antarctic Circumpolar Current (Hausmann & Czaja, 2012;Hausmann et al., 2017).
To date, while THFF is known to increase towards smaller scales, the smallest spatial scale used to quantify THFF is ∼100 km.
The magnitude of THFF depends on the adjustment of the atmospheric boundary layer (ABL) to the SST anomaly.It is suggested that the removal of heat by surface winds is a key process (Bretherton, 1982;Hausmann et al., 2016).On smaller scales, atmospheric heat anomalies are quickly advected away from the SST anomaly, maintaining a large air-sea temperature contrast and strong heat flux damping.While on basin scales, heat advection becomes less efficient (slower), resulting in a small temperature contrast and reduced damping.On global scale, this adjustment completely disappears: the heat removal is controlled by radiation out to space and the THFF reaches only about 1-2 10.1029/2021GL095407 2 of 10 Hayes et al., 1989;Putrasahan et al., 2013;Small et al., 2019;Wallace et al., 1989).A warm mesoscale SST anomaly transfers heat through turbulent heat fluxes up into the ABL.This heat addition reduces stability, enhances vertical mixing, and reinforces the downward transfer of momentum, strengthening surface winds.The opposite occurs over a cold SST anomaly.Past research on mesoscale air-sea exchanges largely focuses on momentum fluxes that is Renault et al. (2016Renault et al. ( , 2019) ) and Seo et al. (2016).However in eddy-rich regions, mesoscale-induced air-sea turbulent heat fluxes play an important role in altering eddy kinetic and potential energy and dampening SST anomalies (Bishop et al., 2020;Ma et al., 2016).Furthermore, mesoscale SST-turbulent heat flux exchanges can strengthen western boundary currents (WBC) by 20%-40% and weaken thermal stratification in the upper ocean (Ma et al., 2016;Shan et al., 2020;Small et al., 2020).It is therefore important to quantify THFF over transient mesoscale eddies.
Observational estimates of THFF at the oceanic mesoscale are restricted by the availability of high-resolution ocean and atmosphere data.First, the consistency and effective resolution of global air-sea heat flux datasets are questionable, due to the different space-time resolutions from either atmospheric reanalysis or satellites (Cronin et al., 2019;Leyba et al., 2016;Li et al., 2017;Tomita et al., 2019;Villas Bôas et al., 2015).Second, the radii of observed mesoscale eddies maybe be overestimated by a factor of 2 due to the interpolation of along-track sea surface height measurements by satellite altimeters into regular grids (Chelton, 2013;Cronin et al., 2019;Ducet et al., 2000;Hausmann & Czaja, 2012;Minobe et al., 2008;Moreton et al., 2020;Small et al., 2008;Xie, 2004).As a result, this study uses a global coupled climate model with higher spatial ocean and atmospheric resolution than currently available in observations.Current state-of-the-art climate models can provide global eddy-rich ocean simulations, with a horizontal resolution of ∼1/ 12 E  .At this resolution, mesoscale eddies can be explicitly resolved globally, except in the highest latitudes with more, smaller and longer-lasting eddies compared to a 1/ 4 E  resolution (Haarsma et al., 2016;Hewitt et al., 2017;Moreton et al., 2020;Roberts et al., 2019).However, whether an eddy-rich ocean results in an improved representation of mesoscale SST-turbulent heat flux exchanges remains to be determined.The ratio of ocean-atmosphere horizontal resolution is likely to be an important factor (Jullien et al., 2020;Wu et al., 2019).In many current high-resolution coupled models with a NEMO ocean component, air-sea fluxes are computed on the atmospheric grid, which requires the interpolation of SST from the oceanic grid to the often coarser atmospheric grid through the OASIS3-MCT coupler (Williams et al., 2018;Valcke, 2013).The interpolation is likely to smooth out mesoscale features resolved on the ocean grid before calculation of the air-sea exchanges and if so, to introduce significant biases in air-sea feedbacks.
Therefore, our study has two goals: (a) to provide the first estimate of THFF over coherent mesoscale eddies globally at smaller spatial scales than previously evaluated and (b) to evaluate if THFF is dependent on the ratio of ocean-atmosphere resolution in coupled models.The estimates are obtained for coupled eddy-resolving and eddy-permitting simulations from the HadGEM3-GC3.1 model.The configurations and methods are introduced in Section 2. Section 3 presents the results addressing the two goals, and Section 4 concludes and discusses implications for future research and model development.

Model Data
We use output from the high-resolution global coupled climate model, HadGEM3-GC3.1 (Williams et al., 2018).The model simulations follow the CMIP6 HighResMIP protocol, as part of PRIMAVERA (Haarsma et al., 2016;Roberts et al., 2019).Three configurations with a different ratio of ocean-atmosphere resolution are compared: N512-12 ( E  25 km atmosphere, 1/ 12 E  ocean), N216-12 ( E  60 km atmosphere, 1/ 12 E  ocean) and N216-025 ( E  60 km atmosphere, 1/ 4 E  ocean).Model outputs are obtained after a 20-year spin-up, and one year of daily data is used (the results are independent of the year chosen).
To compute air-sea latent and sensible heat fluxes, the OASIS3-MCT coupler passes the ocean model SST to the atmospheric grid using a second-order conservative interpolation (Hewitt et al., 2011;Valcke, 2013;Valcke et al., 2015).Here, we define the turbulent heat fluxes (THF) as the sum of latent and sensible heat fluxes, using the convention that positive THF denotes fluxes upwards from the ocean to the atmosphere.In the following, surface air temperature is taken at 1.5 m and the SST on the ocean grid ( O E SST ) is distinguished from the regridded SST on the atmospheric grid ( A E SST ). 10.1029/2021GL095407 3 of 10

Eddy Tracking and Compositing
From SSH outputs from the model simulations, closed coherent mesoscale eddies are identified and tracked daily in the global ocean for 20 years from SSH, using an eddy tracking algorithm adapted from Mason et al. (2014), which is originally based on Chelton et al. (2011).
Briefly, the algorithms detects closed SSH contours around SSH maximum/minimums.Eddy detection is also subject to certain criteria such as a shape test, that is how circular an eddy is.For further details of how the algorithm works, its adaptations and a model comparison to observations, the reader is referred to Moreton et al. (2020).The latter also provides a comparison with altimeter-based results (Ducet et al., 2000).It shows that the observational product likely overestimates the eddy radii because of the processing involved in generating a gridded data set from the satellite tracks.
To isolate mesoscale anomalies, a 10-year climatological mean is removed from the fields, which are subsequently high-pass filtered, by removing a low-pass field obtained by a Gaussian filter of widths 20 E  (zonal) by 10 E  (meridional) (same filter as applied to the SSH for eddy tracking; see Supporting Information S1 for details).Following Frenger et al. (2013), Hausmann and Czaja (2012), Villas Bôas et al. (2015), "composite averaging" is used to remove high-frequency variability associated with weather.Highpass filtered anomalies centered on each eddy are first resized by the effective eddy radius eff E L before averaging.eff E L is defined as the radius of a fitted circle with the same area as the outermost closed SSH contour in each tracked eddy.Rotating the anomalies (to align with background SST or wind direction) before averaging makes little difference to our results.
Finally, the eddies and their associated fields are binned according to their eddy amplitude E A , defined as the absolute difference between either the maximum (anti-cyclones) or minimum (cyclones) SSH and the value of the outermost closed SSH contour of the tracked eddy, from 3 E  0.05 cm (small-amplitude) to 34 E  6 cm (large-amplitude).A global map of the averaged E A per 1 E  squared is shown in the Figure S1 in Supporting Information S1.As expected, larger amplitude eddies are concentrated in eddy-rich regions, such as WBCs and the Southern Ocean.The number of eddy snapshots in each amplitude bin is given in Table S1 in Supporting Information S1.SST and THF from large-amplitude eddies, while a replica for small-amplitude eddies is found in the Figure S2 in Supporting Information S1.Stippling indicates values which are not statistically significant from zero (using student's t-testing with a 99% confidence level).Note that closed contours of the composite anomaly are found beyond one eff E L : this is because eff E L is identified on individual eddies, while the composite averages remove much of the noise revealing close contours beyond eff E L .It is noted that eddy amplitude and eddy radius are not strongly related (Chelton et al., 2011;Moreton et al., 2020).Instead, eddy amplitude ( E A  25 cm) is linearly related to SST anomalies, as shown in Figure S3a in Supporting Information S1 and in previous studies (Villas Bôas et al., 2015).
An accurate comparison of eddy composites from the model to observations is difficult, due to the coarser resolution found in observations and differences in either how the SSH anomalies are isolated (i.e., by standard deviation of SSH anomalies or eddy tracking), the eddy tracking algorithm or the scales retained in the high-pass filtering.Despite this, the O E SST composites in the model have similar magnitudes and spatial distributions to previous observational studies (Frenger et al., 2013;Gaube et al., 2015;Hausmann & Czaja, 2012;Sun et al., 2020).For all resolutions, maximum SST anomalies of E  0.6 K are found in eddies of amplitude of 15 cm (i.e., in eddy-energetic regions, Figure S3a in Supporting Information S1), close to the value of 0.75 K seen in observations (Hausmann & Czaja, 2012).

Decomposition of the Turbulent Heat Flux Feedback
The THFF E  is defined as: where primes indicate the high-pass filtered anomalies, and < > .indicates the eddy-centric composites computed for all eddies tracked in the SSH model outputs.A positive value of E  represents a negative heat flux feedback, that is a dampening of the SST anomaly by the THF.
Due to the regridding of SST to calculate air-sea heat fluxes in the coupled model, two THFFs can be computed from either The THFF O E  relates the THF anomalies to the prognostic SST anomalies in the ocean component, while A E  represents the THFF after re-gridding the ocean grid SST to the atmospheric grid (  does not affect directly the prognostic state of the simulation.
To understand the behavior of the THFFs First, the THF restoring coefficient A E  is a simplification of the latent and sensible heat flux (LHF and SHF) bulk formulae used in the model (Large & Yeager, 2004).Following Frankignoul et al. (1998)   ) or based on re-gridded SST ( A E  ), we can provide an estimate for how the THFF is affected by the ratio of ocean-atmosphere resolution in coupled models.By re-arranging Equations 4-6, relationships between the coefficients can be derived, in order to trace changes from the THF restoring coefficient The THFF  by the air temperature adjustment in the ABL (Equation 7).When the ABL temperature adjustment is weak (i.e.,  is close to the restoring embedded in the THF bulk formulae (i.e., A E  here).Whilst when the adjustment is strong, the THFF A E  , and subsequently the dampen- ing of SST anomalies, is much smaller than predicted by A E  (Frankignoul et al., 1998).In other words, the coefficient A E  represents an upper bound for A E  , which is achieved when air temperature adjustment ( E  ) is zero.This upper bound is the "fast limit" discussed by Hausmann et al. (2017).
The THFF using ocean model SST ( It is anticipated that g E R is smaller than 1 and therefore that In practice, the above coefficients are estimated over coherent mesoscale eddies through linear regressions between data from the composite maps.To remove variability occurring outside the detected eddies (Figure 1), only data within a square of 2 eff E L  2 eff E L is used in the linear regressions.Sensitivity to this choice will be discussed.Regressions for anti-cyclonic and cyclonic eddies are calculated separately, and a weighted average is calculated, using the number of anticyclonic and cyclonic eddies, to produce a total value (given as text in Figure 2).The gradients of linear regression are dependent on SST being on the E x -axis. 10.1029/2021GL095407 5 of 10 Assuming a normal distribution of data and using the student's t-test, 95% confidence intervals are supplied in Figures 2 and 3.

First the THFF coefficients, A E
 and O E  , are discussed for the N512-12 configuration.This configuration is presented first because it is the least affected by regridding biases (Section 3.1).A comparison to N216-12 and N216-025 configurations follows, to evaluate the impact of changes in the ratio of ocean-atmosphere resolutions on the THFF (Section 3.2).

Estimating THFF Over Large-Amplitude Mesoscale Eddies
Figure 2 illustrates the relationships between the composite fields for the large amplitude eddies ( E A = 34 E  6 cm) globally in N512-12. Figure S4 in Supporting Information S1 shows the corresponding plots for small-amplitude mesoscale eddies ( There is a strong linear relationship between the composite anomalies of THF and air-sea temperature contrast (Figure 2a).This supports the linearization of LHF underlying Equation 4 (further supported is provided by a 0.98-0.99correlation between SST and the 1.5 m specific humidity air E Q over coherent eddiesnot shown).The robust estimate of A E  at Frankignoul et al. (1998) and Rahmstorf and Willebrand (1995) and the upper bound of 25-35 Hausmann et al. (2017).This discrepancy could reflect differences in the estimation methods.Published estimates are based on the linearization of bulk formulae using constant drag coefficients and monthly mean large-scale winds.In contrast, our estimates (Figure 2a) implicitly account for (a) the full complexity of the bulk formulae implemented in HadGEM3-GC3.1,where the drag coefficient is function of ABL stability The estimates combines cyclonic and anticyclonic eddies as described in Section 2. In subplots c and e, the regression lines for anticyclonic and cyclonic eddies are plotted in red and blue respectively. 10.1029/2021GL095407 6 of 10 and surface winds (Hewitt et al., 2011) and (b) dynamical adjustments in the ABL such as the modulation of surface winds by mesoscale eddy SST anomalies (Frenger et al., 2013;Roberts et al., 2016).
The atmospheric adjustment parameter E  is estimated at 0.34 E  0.01 for large amplitude eddies (Figure 2b), that is the surface air temperature air E T anomaly is about a third of the mesoscale SST anomaly.Previous studies give 0.5 in the WBCs and the Antarctic Circumpolar Current (ACC) core, increasing to 0.9 in quiescent regions (Hausmann et al., 2017).However, these estimates are limited by the scale of ERA-I reanalysis (0.75 0.75 E   ) and do not isolate coherent eddies.Although the modeled large-amplitude eddies used in Figure 2 are mostly found in WBCs (Figure S1 in Supporting Information S1), our estimate suggests that air E T adjustments drop further below 0.5 over coherent mesoscale eddies.2d) can now be explained by combining estimates of A E  and E  using Equation 7:

The value of
As most large-amplitude eddies are found in the WBCs, our modeled estimate of A E  agrees well with previous observational estimates of 40-56 Kuroshio region and 40 (Hausmann et al., 2016;Ma et al., 2015).Finally, the THFF on the prognostic SST, O E  , is about 25% smaller than 2e).The reduction reflects the 25% decrease in the amplitude of mesoscale SST anomalies brought by the SST regridding ( 0.74 , see Equation 8; Figure 2c). is significant, and can be at- tributed to the regridding between 2).This results in an uncertainty in W m K E   (found consistently across all eddy amplitudes, and all resolutions).Interestingly, a small asymmetry between cyclonic and anticyclonic eddies in O E  can also be attributed to g E R (Figure 2 and Figure S4 in Supporting Information S1), potentially due to slight differences in magnitude of the eddy anomaly.It therefore appears that the regridding, even in the most favorable case of near matching resolutions, is a source of noise and non-linearities.Figure 2 is repeated in Figure S5 in Supporting Information S1 using data from the whole composited region shown in Figure 1, that is a 5.6 5.6 square.The asymmetry between polarities vanishes, which suggests this is not a robust feature, but possibly an artifact from the tracking algorithm and/or the regridding process.We do not investigate this asymmetry further. developed above for large-amplitude eddies applies equally well to small-amplitude eddies (see Figure S4 in Supporting Information S1).We therefore present variations of

The rationalization of the THFF
 as a function of eddy amplitude E A in N512-12 (Figure 3a).To first order, the THFF increases with eddy amplitude (and hence with mesoscale SST anomalies, see Figure S6 in Supporting Information S1).From a minimum THFF of Referring to Equation 7, variations in  roughly increases with the eddy amplitude, or equally with the eddy SST anomaly as the two are strongly correlated (see Figure S3a in Supporting Information S1).This likely reflects non-linearities embedded in the bulk formulae.One such non-linearity is the effect of the surface wind speed.As highlighted in previous studies (e.g., Roberts et al., 2016, and references therein), the ABL response to mesoscale SST anomalies includes a surface wind speed response proportional to the mesoscale SST anomalies.Here, we confirm that, as expected, the wind speed anomaly increases with the eddy amplitude (Figure S3b in Supporting Information S1).This effect contributes to strengthen the air-sea exchanges A E  over large eddies.However, it is likely that other non-linearities play a role (as suggested by results for other configurations, see below). )with eddy amplitude is slightly smaller, compared to N512-12.This is consistent with a weaker surface wind response in N216-12 and N216-025 (Figure S3b in Supporting Information S1).The near absence of a surface wind response in N216-025 suggests that other non-linearities such as the dependence of drag coefficient on temperature and ABL stability, contribute to the dependence of  exhibits a strong linear correlation with the regridding parameter g E R (Figure 3d), with a slope of E  1 as predicted by our simplified relationships (see Equation 8).This reinforces our interpretation that the regridding of SST (captured by

Conclusions
Turbulent heat flux feedback over coherent mesoscale eddies is estimated globally in three configurations of the high-resolution coupled model HadGEM3-GC3.1.First, for the highest ocean-atmosphere resolution available (where the impact of SST regridding from the ocean grid to the atmosphere grid is minimal), the estimates of the THFF over mesoscale eddies range from 35 to 45 where values roughly increase with eddy amplitude.This is the first time such estimate is provided as previous studies did not resolve such small scales nor attempted to isolate coherent eddies.Second, we investigate configurations with larger mismatch between oceanic and atmospheric resolutions.We find that the regridding of SST from the ocean to atmosphere grid can underestimate the eddy-induced THFF by 20%-80%.Importantly, this low bias increases with the ratio between atmospheric and ocean resolutions, implying that increasing the oceanic resolution at constant atmospheric resolution can actually degrade solution, at least in the representation of air-sea feedbacks.
The low bias in the THFF suggests that eddies are not dampened enough in the model.Eddies have a first order impact on the dynamics of WBCs and the ACC.However, small-amplitude eddies that dominate the eddy population cover the global open ocean, influencing the stratification, ocean heat uptake and biological processes.These eddies have a strong THFF of 35-40 and are the most affected by the low biases due to regridding.Further work is needed to understand these biases, but it is likely to have range of impacts beyond eddy-rich regions: artificially large SST anomalies are likely to cause an artificially large local and large-scale ocean and atmospheric response (Bishop et al., 2020;Frenger et al., 2013;Ma et al., 2016).
Our findings should be tested with other high-resolution climate models, which adopted different coupling strategies (Yang et al., 2018).In addition, while our focus was on horizontal resolution, it is likely that the vertical resolution, in both the ocean and atmosphere, play a major role in the representation of mesoscale air-sea exchanges through its influence of the ABL adjustment (Stewart et al., 2017).Finally, we leave binning by eddy radii and exploring the effect of lags between SST and THF on our THFF estimates for future work.
The results in this study hold implications for future model development.Similarly to HadGEM3-GC3.1,many current high-resolution coupled models (which use the OASIS coupler for example) compute air-sea turbulent heat fluxes on the atmospheric grid, using regridded SST (Roberts et al., 2019;Valcke et al., 2015).For the long spin-ups needed for climate simulations, it is unrealistic to expect the atmospheric resolution to match the oceanic resolution.Instead, it is advised when resolving mesoscale eddies, that air-sea heat fluxes should be calculated on the finer-scale oceanic grid, as done by the Community Earth System Model (see Yang et al., 2018).This method ensures that the high-resolution SST anomalies are maintained, although this requires a large logistical change for many coupled models and is computationally much more expensive.Our results also indicate that the regridding introduces a noise and an asymmetry between cyclonic and anticyclonic eddies.Essentially, we need a "better" regridding of O E SST to A E SST although it is inevitable that even the best regridding technique will degrade mesoscale SST anomalies in large ocean-atmosphere resolution difference.In ocean-only models, the ocean component is driven through bulk formulae and prescribed surface atmospheric fields, that is without ABL adjustment (i.e., 0 E   in our notations).In such setups, we expect mesoscale THFF to approach A E  .However, the absence of an ABL adjustment also influ- ences A E  (e.g., neglecting the effect of dynamical adjustment on the drag coefficient).The net effect of these assumptions on the mesoscale THFF in ocean-only models remains to be quantified.

Figure
Figure1shows composites of

Figure 1 .
Figure 1.Composite maps of turbulent heat flux (THF) in W m K   2 1


and Hausmann et al. (2017), we assume that the LHF can be linearized to be expressed in terms of the air-sea temperature difference, equals zero there is no ABL response or adjustment, whilst when E  equals one, a complete adjustment occurs resulting in a zero THF.Third, the g E R coefficient measures the impact of the ocean-to-atmosphere regridding on the SST magnitude.If g E R equals one, the magnitude of the SST anomalies is preserved during the regridding.By isolating THFF based onO E SST ( O E

Figure 2 .
Figure 2. Relationships between the composite fields of SST O A / , turbulent heat flux and air E T , with the estimated coefficients ( / O A


are mainly driven by changes in the THF restoring A E  whilst the atmospheric adjustment E  is relatively insensitive to eddy amplitude (compare Figures S3d and S3e in Supporting Information S1).The restoring coefficient A E

Figure
Figure 3 summarizes estimates ofA


depends greatly on the difference between the oceanic and atmospheric grid resolutions: -40% for N512-12, to 40%-60% for N216-025 and to approximately 60%-80% for N216-12.Crucially, the low bias is the largest for the smaller amplitude eddies, which cover most of the global ocean in the configuration with the largest ratio between atmospheric and oceanic resolutions, N216-12.The typical eddy scale of small amplitude eddies ( eff E L  40 km on average) is smaller than the 10