Impact of resolution on the atmosphere–ocean coupling along the Gulf Stream in global high resolution models

We have investigated the horizontal resolution dependence of the ocean–atmosphere coupling along the Gulf Stream, of simulations made by six Global Climate Models according to the HighResMIP protocol, and compared it with reanalysis and remote sensing observations. Two ocean–atmosphere interaction mechanisms are explored in detail: The Vertical Mixing Mechanism (VMM) associated with the intensification of downward momentum transfer, and the Pressure Adjustment Mechanism (PAM) associated with secondary circulations driven by pressure gradients. Both VMM and PAM are found to be active even in the eddy-parameterized models. However, increasing ocean and/or atmosphere resolution leads to enhanced ocean–atmosphere coupling and improved agreement with reanalysis and observations. Our results indicate that while one part of the stronger air–sea coupling is attributable to the refinement of the oceanic component to eddy-permitting, optimal results are obtained only by further increase of the atmosphere resolution too. The use of the eddy-resolving model show weaker or same coupling strength over the eddy-permitting resolution. We conclude that at least eddy-permiting ocean resolution and comparable atmosphere resolution are required for a reliable ocean–atmosphere coupling along the Gulf Stream.


Introduction
There is growing evidence that strong sea surface temperature (SST) gradients, oceanic fronts and eddies can force the atmosphere by affecting the near surface wind. At basin scales there is a negative correlation between SST and surface wind speed. Higher wind speeds are located over cooler SSTs, as the atmosphere is forcing the ocean through evaporative cooling and entrainment of colder thermocline waters to the surface (Cayan 1992;Liu et al. 1994;Alexander et al. 2002;Deser et al. 2010). In contrast, at smaller scales, in regions of large mesoscale oceanic activity, higher wind speeds are associated with warm SST anomalies O'Neill et al. 2005;Xie 2004). This reversed relation indicates that at mesoscale, the ocean is forcing the atmosphere. The surface wind response to the spatial variability of SST generates divergence and curl of the near surface wind (Chelton et al. 2004). The proposed potential mechanisms, which are responsible for the SST influence on surface winds, are the VMM and the PAM. Additionally, there are recent studies suggesting that atmospheric fronts 1 3 and storm systems can also affect and contribute to the timemean near-surface wind convergence (Parfitt and Seo 2018;O'Neill et al. 2017;Parfitt and Czaja 2016), but this contribution is beyond the scope of this study.
The VMM describes that over warm SST anomalies, there is an intensification of the downward momentum transfer due to the destabilization of the overlying marine atmospheric boundary layer (MABL), which results in acceleration of surface winds (Wallace et al. 1989;Sweet et al. 1981). Acceleration of the wind over warm waters, where it is blowing parallel to the SST front, results in horizontal wind shear, leading to low-level vorticity. In contrast, acceleration where the wind is blowing across the SST front, results in wind stress divergence (convergence) when the wind flow is from a cold to a warm (warm to cold) area (Chelton and Xie 2010;Chelton et al. 2004;O'Neill et al. 2003). According to this, the downwind (crosswind) SST gradients are proportional to the wind stress divergence (curl) (Chelton et al. 2001(Chelton et al. , 2004O'Neill et al. 2003).
On the other hand, PAM describes the response of secondary circulations in the MABL which are driven by pressure gradients due to the differential heating of the atmosphere on each side of the front, leading in turn to surface wind stress convergence (divergence) over warm (cold) SST patterns (Hsu 1984;Lindzen and Nigam 1987;Warner et al. 1990;Wai and Stage 1989). We can visualize this mechanism by the analogy of the sea-breeze circulation proposed by Hsu (1984), according to which the warm ocean along the Gulf Stream (GS) has the role of the land. Minobe et al. (2008) noticed that the wind divergence pattern co-locates with the Laplacian of sea-level pressure ( ∇ 2 SLP), and that both have many similarities with the pattern of the signreversed Laplacian of SST ( −∇ 2 SST).
The importance of these mechanisms is not limited to the fundamental understanding of the Earth system. The SST is not only modifying the near surface wind fields but also can affect the free troposphere influencing that way the Northern Hemisphere from weather to climate time-scales (Minobe et al. , 2010Ma et al. 2015;Kobashi et al. 2008;Lee et al. 2018, Roberts et al. in press). This is through the wind convergence in the MABL that anchors a band of precipitation and deep upward motion extending to the free troposphere, with a remote impact by forcing of atmospheric planetary waves . Moreover, SST induced modifications of the wind speeds can affect cyclogenesis, and influence the storm tracks and storm intensity (Xie et al. 2002;Nakamura et al. 2004;Small et al. 2014;Woollings et al. 2010). In addition, the near surface wind stress response to the SSTs can possibly feed back on the ocean by altering, for example, Ekman pumping (Chelton et al. 2007;Spall 2007).
Due to the coarse resolution of climate models, we hypothesize that the impact of the ocean on the atmosphere at small scales is not correctly represented, with implications for the atmospheric dynamics such as position and strength of the storm tracks (Piazza et al. 2016;Parfitt et al. 2017;Small et al. 2014Small et al. , 2019 and that higher resolution is needed. Recently, a new set of global high resolution climate simulations following the HighResMIP protocol (Haarsma et al. 2016) have become available. Using these simulations, we investigate the ocean-atmosphere coupling in the GS region that is characterized by sharp fronts and large eddy activity. This region has also been the focus of earlier studies of VMM and PAM (e.g. Wai and Stage 1989;Minobe et al. 2008;Chelton and Xie 2010). Air-sea interaction over the GS region is of great interest since it affects cyclogenisis and storm track characteristics (Small et al. 2014), induces deep tropospheric response and influences the large scale atmospheric circulation . The High-ResMIP protocol, requiring all simulations to be performed with a low and high resolution model configuration, allows us to investigate the spatial resolution dependence of the ocean-atmosphere interaction. The impact of increased resolution in HighResMIP simulations on climate and variability has been documented in several studies, including Moreno-Chamarro et al. (2021) who discuss the impact on biases. We have focused here on the VMM and PAM mechanisms and compared the results with available observations. In the following Sect. 2 we describe the models' configurations, the observation/reanalysis datasets and the techniques used for this analysis. In Sect. 3 we present our results, followed by a discussion in Sect 4. Finally, a summary of the conclusions is given in Sect. 5.

Data
In the current study, we analyse the "control-1950" simulations of the High Resolution Model Intercomparison Project (HighResMIP) (Haarsma et al. 2016), which is part of phase 6 of the Coupled Model Intercomparison Project (CMIP6) (Eyring et al. 2016). This experiment consists of 100-year coupled simulations of different horizontal resolutions, with external forcings kept constant at 1950's values. All the simulations are forced with greenhouse gases (GHG) and aerosols loadings using 1950s climatology (fixed values). The initialization of the coupled runs was provided by ∼50-year spin-up integrations starting from the observed ocean and atmosphere state in 1950, with fixed forcing representative of the 1950s. Details of the experimental set-up of the HighResMIP spin-up and control-1950 simulations are described in Haarsma et al. (2016).
The six GCMs examined in this study are from the EU funded Horizon-2020 PRIMAVERA project. Each of them consists of more than one configuration (15 in total). The differences between the configurations of each model are appertain to the horizontal resolution of the atmosphere, ocean or both components. Resolution-dependent parameterizations might differ among the configurations. Particularly, the Gent and McWilliams (1990) scheme for the ocean eddy parameterization is deactivated in the HadGEM3-GC3.1, ECMWF-IFS and CMCC-CM2 configurations of oceanic resolution 25km or higher. For each model configuration, one member was analyzed. We argue that for the ocean-atmosphere coupling along the GS, natural variability is sufficiently sampled in a single 100-year integration. This is supported by the fact that applying our analysis to only the first 50 years and only the last 50 years of the simulations, yields similar results (Supplemental Fig. S1). The model configurations and their characteristics are summarized in Table 1. Detailed descriptions of the configurations can be found in Roberts et al. (2019) for HadGEM3-GC3.1, in Roberts et al. (2018b) for ECMWF-IFS, in Haarsma et al. (2020) for EC-Earth3P, in Voldoire et al. (2019) for CNRM-CM6, in Gutjahr et al. (2019) for MPI-ESM1-2 and in Cherchi et al. (2019) for CMCC-CM2. A short description of these models is provided in the Supplemental Material (text S1).
For the atmospheric components, the effective resolution in addition to the nominal resolution is listed in Table 1 as computed in Klaver et al. (2020). The effective resolution is based on the shape of the kinetic energy spectrum and provides an estimate of the boundary of the simulated scales that are not affected by the resolution of the model. Satellite observations combined with reanalysis data, for the period January 2007 to December 2018, are used for comparison with the model output. This period is chosen because of the availability of high resolution satellite observations. The 11-year observational period is, however, much smaller than the 100-year control simulation which can increase the uncertainties. Below we briefly describe the three different observational and reanalysis data-sets used in this study.
SST measurements are obtained from the National Oceanic and Atmospheric Administration (NOAA) 0.25 • daily Optimum Interpolation SST version 2 (OISST.v2) (Banzon et al. 2016). This data-set is constructed by combining observations from the Advanced Very High Resolution Radiometer (AVHRR) and other platforms (ships, buoys) on a regular global grid. Optimum interpolation (OI) is used to produce a spatially complete map.
Wind stress divergence and curl are obtained from Advanced SCATterometer (ASCAT) on MetOp-A (Figa-Saldaña et al. 2002). It is a product of the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Ocean and Sea Ice Satellite Application Facility (OSI SAF) provided through the Royal Netherlands Meteorological Institute (KNMI). We use daily values from the L3 REP product, sampled on a 0.25 • grid.
SLP is derived from the European Centre for Medium Range Weather Forecasts Reanalysis 5 atmospheric ReAnalysis (ERA5) (Hersbach et al. 2020). Data processing for ERA5 is carried out by using ECMWF's Earth System

Methods
We analyze the three-month boreal winter season from December to February (DJF). The rationale for this choice is that due to the increased temperature difference between the atmosphere and ocean, and the existence of stronger winds, the air-sea coupling is stronger in winter (Minobe et al. 2010;Putrasahan et al. 2013).
We use monthly values of the surface wind stress components, SST and SLP from the 100-year control simulations. We choose to work with monthly values in order to remove the effects of energetic synoptic weather variability (Song et al. 2009;Chelton et al. 2007;Chelton and Xie 2010). Since we are investigating two mechanisms which explain the atmospheric response to the ocean, we first bilinearly interpolate the SST onto the atmospheric grid of each model configuration. In that way, our computations are not affected by the oceanic variability at spatial scales that cannot be resolved by the atmosphere. Following O'Neill et al. (2003), we compute the downwind SST gradient as the projection of the SST gradient on the wind stress vector. The crosswind SST gradient is the component of the SST gradient that is perpendicular to the projection of the SST gradient on the wind stress vector. In addition, the wind stress divergence and curl, the ∇ 2 SLP and the −∇ 2 SST are computed. For the calculation of the derivatives, a central difference scheme is used.
We compute the detrended anomalies of the fields. Next, following Roberts et al. (2016), we applied a high-pass spatial filter to these anomalies, by subtracting smoothed fields that are computed by the application of a 18 • longitude × 6 • latitude box car filter from the original fields. This filter is similar to a loess spatial high-pass filter with half-power points at 30 • longitude and 10 • latitude which was applied in previous studies to isolate the mesoscale signal (Bryan et al. 2010;Maloney and Chelton 2006).
The linear relationship, in terms of least squares regression, is investigated between the anomalies of the three aforementioned pairs of fields ([wind stress divergence, downwind SST gradient], [wind stress curl, crosswind SST gradient], [ ∇ 2 SLP, −∇ 2 SST]). The slope of the linear fit is the so-called coupling coefficient, as its magnitude indicates the strength of the air-sea coupling. To compute the coupling coefficients, we create monthly binned histograms of the SST-related fields in the GS domain (30 • -50 • N, 80 • -40 • W, boxed domain in Fig. 1). Then, we create their monthly histograms by determining the averaged atmospheric fields for each bin of the associated SST-related fields. Next the monthly histograms over the analyzed time period are averaged and plotted in a scatter diagram. The linear fit is made over the mean histogram, resulting in one coupling coefficient for each configuration.
For a fair comparison with the observations, we handle the latter in the same way as the model data. The only difference is that for the ASCAT products we use daily values of the wind stress divergence/curl and then we compute the monthly averages for the points with 20 or more available observations during a month. We proceed this way, since due to the inhomogeneous distribution of the observations over space and time, the fields appear noisy and this 'noisy' character is further amplified by the derivatives operation. The coupling coefficient though, is not affected significantly by first averaging the winds stress and then computing the divergence (Supplemental Fig. S2).
The use of the monthly mean values in our analysis could be considered a limitation. As we do not compute the variables of interest with higher time-resolved data (e.g. daily) before averaging, errors may occur in the non-linear terms of the downwind and crosswind SST gradients. For this reason, we also compute the coupling coefficient by first applying the operators on the daily values before averaging, for the configurations of the HadGEM3-GC3.1 model and the observations, in order to investigate if any significant differences are induced. Any regridding, required before the averaging, was made using a conservative interpolation method. Even though the magnitude of the coupling differs between the two analyses, the dependence of the coupling on the horizontal resolution is not sensitive to the order of the aforementioned steps (not shown here).

Results
The 11-year winter mean climatology of ASCAT near surface wind stress divergence and NOAA OISSTv.2 SSTs are shown in Fig. 1. As previous studies have shown (Minobe Chelton et al. 2004;Takatama et al. 2012), observations reveal wind stress divergence (convergence) along the cold (warm) flank of the GS. Strong divergence (red in Fig. 1) is observed in the region of Grand Banks where the GS turns to the north. Westerlies cross the GS resulting in divergence due to the wind speed differences on each side of the strong SST gradients.
The ability of the models to reproduce these wind features and their dependence on spatial resolution is shown in Fig. 2, where the 100-year DJF climatology of wind stress divergence and SST contours are presented for the different configurations of the models. Generally, all models are capable of capturing the divergence and convergence zones on the cold and warm flank of the GS respectively. However, configurations with a low-resolution oceanic component of about 1 • , where the ocean eddy dynamics are parameterized, display much lower absolute values of divergence. Especially east of Newfoundland, the divergence values are smaller by a factor of four or more compared to the observations. This is associated with the much weaker SST gradients of the eddy-parameterized configurations (Supplemental Fig. S3). By increasing the oceanic resolution from 1 • to 0.4 • -0.25 • , representing eddy-permitting dynamics, the strong SST gradients at about 50 • W around the Grand Banks are better resolved, in line with Kirtman et al. (2012). However, a substantial underestimation of the divergence field is observed even in eddy-permitting models. Even though there is a large spread among the models, enhancing the atmospheric resolution improves the simulated divergence field. The HadGEM3-CG3.1-HH, HM, MM and CMCC-CM2-VHR4 are the configurations with the best performance, with divergence patterns and magnitudes comparable to observations.
In observations, the GS separates from the US coast near Cape Hatteras. Only HadGEM3-CG3.1-HH with an ocean eddy-resolving resolution of 1/12 • and atmosphere resolution of 50 km correctly simulates this separation. Coarser resolutions (oceanic and atmospheric) result in separations further to the north along the eastern coast of North America. This is a well-known feature and has left: eddy-parameterized models, eddy-permitting models, eddy-permitting models with further increase of the atmospheric resolution, eddy-resolving models been addressed in more detail in previous studies (Bryan et al. 2007;Chassignet and Xu 2017;Seidov et al. 2019). Consequently, the associated convergence zones of wind stress are also shifted northwards for coarser resolutions.

Spatial anomalies patterns and correlation
To illustrate the resolution dependence of downwind SST gradient and near surface wind stress divergence and their collocation, we present in Fig. 3 (left) the high-pass filtered downwind SST gradient and the near surface wind stress divergence anomalies for the different configurations of the HadGEM3-GC3.1 model and observations (ASCAT-OISST). For all configurations, the relationship between the two fields is evident as wind stress divergence anomalies are associated with positive downwind SST gradient anomalies, while convergence anomalies are found over negative downwind SST gradient anomalies. The spatial covariability of the fields is large for all configurations. However, the ocean eddy-parameterized (1 • ) HadGEM3-GC3.1-LL results in lower absolute values of anomalies for both fields by a factor of five, compared to the eddy-permitting (0.25 • ) and eddy-resolving (1/12 • ) configurations. In addition, the eddy-parameterized Means are over one representative December-February period as the simulations do not represent actual years. The contours intervals (c) differ for each resolution for better visualization and are noted (in Pa/10000km) at the left top corner in each case. The colorbar range differs for the LL configuration configuration exhibits patterns of a much larger scale. Configurations with eddy-permitting oceanic resolution, result in increased magnitude of the downwind SST gradient anomalies and smaller spatial features. The collocation of the fields becomes more evident, by further increasing the atmospheric resolution, from 100 to 50 km. Further increase of the oceanic resolution to eddy-resolving, leads to enhanced values of downwind SST gradient anomalies without a substantial associated increase in the wind stress divergence. Similar dependence on resolution is observed for the crosswind SST gradient and wind stress curl covariability as illustrated in Fig. 3 (right).
The other models reveal a similar resolution dependence of the SST influence on near surface winds as the HadGEM3-GC3.1 model (Supplemental Fig. S4-8). All the eddy-parameterized configurations are incapable of reproducing the small scale features. Due to their coarse resolution, it is not possible to resolve the mesoscale activity of the GS resulting in small anomalies of the aforementioned fields. Increasing oceanic and atmospheric resolution results in higher amplitude of anomalies with smaller scale structure and better agreement with observations.
To further investigate the relation between the near surface winds and the SSTs, and its variability in space in the GS region, Fig. 4 shows the temporal correlation between the downwind SST gradient and the wind stress divergence for all models and observations/reanalysis. All models reveal strongest correlation along the strong SST gradient of the GS, as expected. Models with eddyparameterized ocean components (1 • ) exhibit low correlations over a large area. The largest correlations with values larger than 0.5 are obtained at the end of the GS region east of Nova Scotia. Approximately 50% of the values of the correlation coefficient in eddy-parameterized models are lower than 0.2.
The eddy-permitting (0.4 • -0.25 • ) models generally reveal higher correlations up to 0.8 over a much larger area than the eddy-parameterized models. The overall increase of the correlation in the ECMWF-IFS-MR, reveals that the oceanic resolution has a primary role in this improvement. However, Fig. 4 Correlation of timeseries of spatially high-pass filtered wind stress divergence and downwind SST gradient anomalies. Correlations are computed using monthly values for the winter period (DJF).
The stippling shows statistically non-significant correlations at 95% significance level using two-sided student's t-test 1 3 even though there is spread between the models, the correlation appears quite sensitive to the atmospheric resolution too. We note that the band of positive correlation along the GS is further enhanced, when only the atmospheric resolution increases (HadGEM3-HM, CMCC-VHR4, MPIESM-XR and ECMWF-IFS-HR), displaying a 50-100% increase of the area with correlations larger than 0.6. The only eddyresolving (1/12 • ) model HadGEM3-CG31-HH displays lower correlations than the eddy-permitting version with the same atmospheric resolution HadGEM3-CG31-HM. Possible reasons for this unexpected behavior will be discussed below.
In general, ASCAT data display lower correlations in comparison with the models, except for the eddy-parameterized ones. Similar behavior was noticed by Roberts et al. (2016) where the correlation between the wind stress and the SSTs was examined. This can possibly be explained by the temporal and spatial inhomogeneity of the satellite observations. In Fig. 1, we can notice the "noisy" character of the observed divergence field, which is indicative of a certain degree of uncertainty in the observations due to the limited period covered by the data. Figure 5 displays the correlation coefficient between the high-pass filtered crosswind SST gradient anomalies and the associated near surface wind stress curl (note that reversed colors are used compared to Fig. 4). The dependence of the amplitude and spatial structure of correlation on the spatial resolution of the model components is similar to the one found for the downwind component of SST gradient and the divergence of the wind stress. Negative (positive) crosswind SST gradient anomalies are associated with cyclonic (anticyclonic) wind stress anomalies resulting in anti-correlation. The correlation is systematically lower for vorticity than for divergence. The possible reasons for this will be discussed below. In addition, while the HadGEM3-HH shows lower correlations than the HadGEM3-HM for the relationship between the downwind SST gradient and the wind stress divergence, here the two configurations display similar results.

Fig. 5
Correlation of timeseries of spatially high-pass filtered wind stress curl and crosswind SST gradient anomalies. Correlations are computed using monthly values for the winter period (DJF). The stip-pling shows statistically non-significant correlations at 95% significance level using two-sided student's t-test

Coupling coefficient
To quantify the strength of air-sea coupling in the GS region, we focus on the red box area defined in Fig. 1 (30 • -50 • N, 80 • -40 • W). Figure 6 shows the binned scatter plots of spatially high-pass filtered downwind SST gradients and wind stress divergence for the HadGEM3-GC3.1 configurations. The strength of the VMM is determined by the slope of the linear regression. Slopes are defined to be statistically significant when the p-value <= 0.05 , based on a Wald test with t-distribution. The HadGEM3-GC3.1-LL results in a four times weaker slope than observed, s = 0.18 × 10 −2 Pa/ • C. The other three configurations of the HadGEM3-GC3.1 model exhibit much higher coupling with s = 0.52 × 10 −2 Pa/ • C, s = 0.7 × 10 −2 Pa/ • C and s = 0.6 × 10 −2 Pa/ • C for the MM, HM and HH respectively. The lower coupling of HH configuration in comparison with the HM is in agreement with the results in Fig. 4.
The results for all the models and the observations are synthesized in Table 2, displaying not only the coupling coefficients between the wind stress divergence and the downwind SST gradient but also those between the wind stress curl and the crosswind SST gradient. In addition, the the 95% confidence intervals of the slopes are presented, together with the correlation coefficients in parenthesis. All the values are statistically significant.
All the eddy-parameterized models (HadGEM3-LL, ECMWF-LR, EC-Earth-LR and CNRM-LR) exhibit linear relationships with coupling coefficients not exceeding 0.3 ×10 −2 Pa∕ • C between the two pairs of fields. Models of eddy-permitting resolutions of 0.4 • and 0.25 • , give coupling coefficients of larger magnitude. One part of this increase is attributable to the better representation of the mesoscale oceanic activity by the finer oceanic component. This is evident by the comparison of the LR and MR configurations of the ECMWF-IFS, where there is an increase of the coupling coefficient by only refining the oceanic resolution to eddypermitting. However, the CMCC-CM2 HR4 produce coupling coefficients that are similar to the eddy-parameterized models especially in the case of the relationship between the downwind SST gradient and the wind stress divergence. On the other hand, the VMM is sensitive to the atmospheric resolution too. Increasing the atmospheric resolution results in enhanced coupling and leads to values similar in amplitude with the observed for the case of wind stress divergence. The HadGEM3-HH shows weaker or same coupling strength compared to the HadGEM3-HM.
The coupling coefficients between the wind stress curl and the crosswind SST gradient are for all models and their configurations lower than the ones between the downwind component of SST gradient and the divergence of the wind stress. This is in agreement with the observations, but the differences in the models are larger. Reduction is also noticed in the temporal correlations as discussed above. This result is consistent with other studies and has been attributed to many causes. Chelton et al. (2001) and Kilpatrick et al. (2014) attribute this to differences in the Marine Atmospheric Boundary Layer (MABL) response for cross-front and along-front winds. The finite timescale for adjustment of the MABL, allows the boundary layer to come into equilibrium with the SST for an along-front flow, in contrast to a cross-front flow where the downwind SST gradients are sharper (Chelton et al. 2001). That may enhance the wind stress divergence response on the downwind SST gradient  Fig. 1). The calculations were made by averaging the monthly binned histograms for Dec-Feb, over the years 1950-2050. The bin range is 0.1 • C/100km. The error bars represent the ±1 standard deviation over the averaged samples in each bin. The red line is the linear fit in terms of least squares. Statistically significant slopes are represented with green color text. The criterion for a significant value is p ≤ 0.05 (based on a Wald test using t-distribution) due to the increased turbulent fluxes in the non-equilibrium state. Differences in the ocean current speed effects on the cross-isotherm and along-isotherm wind stress can also cause differences in the divergence and curl responses (e.g. Chelton et al. 2004). O'Neill et al. (2010) explain the weaker wind curl response to the SST gradients, by separating the effect of the SST-induced wind speed and direction gradients on the wind curl and divergence. They found that while the SST-induced wind speed anomalies contribute equally to the curl and divergence, the SST-forced wind-direction gradients diminish the wind curl response to crosswind SST gradients and enhances the wind divergence response to downwind SST gradients. Figure 7 shows the observed winter-averaged ∇ 2 SLP (top) and the −∇ 2 SST (bottom), by the ERA5 and the NOAA-OISST products. Negative pressure anomalies are located at the warm flank of the GS while positive pressure anomalies are found over the cold flank. ∇ 2 SLP exhibits many similarities with the −∇ 2 SST and the wind stress divergence (Fig. 1). These results are consistent with Minobe et al. (2008) and suggest the existence of PAM over the GS region, according to which warm (cold) anomalies induce low (high) pressure anomalies and low level convergence (divergence). This is also evident in the models (Supplemental Fig. S9-10). However, the eddy-parameterized models exhibit much weaker −∇ 2 SST and a north-shifted band of positive ∇ 2 SLP, due to the unrealistic separation of the GS. Higher resolution configurations are in better aggrement with observations.

Correlation
Spatial patterns of the temporal correlation between the high-pass filtered −∇ 2 SST and the associated ∇ 2 SLP (Fig. 8), highlight a similar dependence on resolution for PAM as for VMM. The eddy-parameterized models only display correlations larger than 0.5 east of Nova Scotia, due to the weaker SST gradients and the shifted, more zonally oriented path of the GS. The eddy-permitting models show a systematic increase of the correlations, exceeding the observed, with values above 0.5 starting from Cape Hatteras and continue all along the SST front. Enhancing atmosphere resolution leads to higher correlations and improves further the spatial structure. The eddy-resolving HadGEM3-CG3.1-HH shows no improvement compared to the eddy-permitting HadGEM3-CG3.1-HM. A second band of increased correlation is noticed south of 30 • in the extratropics, both in observations and some models. However, the climatological −∇ 2 SST and ∇ 2 SLP do not reveal a signal in that region ASCAT and NOAA-OISSTs 0.75 ± 0.023 (1) − 0.67 ± 0.022 (− 1) (Fig. 7), suggesting that this band is not associated with the PAM mechanism over the GS. We speculate that this band of correlation is a reflection of the enhanced ocean-atmosphere coupling in the tropics (Kushnir et al. 2002).

Coupling coefficient
Since the Laplacian operator acts as a high-pass filter, it is sensitive to the grid on which it is computed. We argue that the linear relationship between the Laplacians is scaledependent, and for this reason a fair comparison between the different resolutions is not possible. In order to demonstrate this we compute the −∇ 2 SST and the ∇ 2 SLP both from ERA5, with the same numerical scheme, first using every grid point (effective resolution of 0.25 degrees, "HR") and subsequently using every fourth grid point (effective resolution of 1 degree, "LR"). The coupling coefficient for the HR is 0.26 while for the LR is 0.41 (Supplemental Fig. S11). This increase of the coupling coefficient, by computing the Laplacians over larger distances using the exact same dataset, indicates that comparing the coupling coefficients between the Laplacians across different model resolutions, is not simply a comparison between the realism of the models. On the other hand, the temporal correlations are not affected by the computation of the Laplacians on different grid spacings (Supplemental Fig. S12). For this reason, we do not proceed to a comparison of the coupling coefficients for the PAM mechanism and we base our conclusions on the temporal correlations, which represent a more reliable metric for this case.

Discussion
Generally, we see that horizontal model resolution plays an important role in the accurate representation of air-sea interactions. Even though, the eddy-parameterized configurations are able to capture the oceanic forcing, they underestimate the coupling strength and exhibit lower temporal/ spatial correlation and coupling coefficients. Increasing the resolution leads to simulations that are in better agreement with the observations. The oceanic resolution appears to be a key feature for more realistic simulations, which is in agreement with other studies (e.g.  2017)). However, only after a further increase of the atmospheric resolution, we obtain results comparable to the observations. Bellucci et al. (2021) investigate the air-sea interaction over the GS region using the same simulations as here. Looking at the correlation and co-variance patterns between the SST and the turbulent heat fluxes (THF) they also find that enhancing the spatial resolution leads to better representation of the coupled ocean-atmosphere processes. However, they attribute the improvements on the finer oceanic resolution, since their results show less sensitivity to the resolution of the atmospheric component. Our results, on the other hand, indicate that for an optimal coupling, in addition to an eddy-permitting oceanic component, the atmospheric resolution should not be much coarser than the ocean resolution. The results from the eddy-permitting models for which are available two atmospheric resolutions (HadGEM3-GC3.1, ECMWF-IFS, MPIESM-1-2, CMCC-CM2), support this hypothesis. Experiments with eddy-parameterized oceanic components coupled with high-resolution atmospheric models could provide a better understanding of the relative importance of the atmospheric and oceanic resolutions, yet such simulations are not available in this set of experiments.
A peculiar result is that for the correlations as well as the coupling coefficients, the eddy-resolving HadGEM3-CG3.1-HH shows lower or similar values compared to the eddypermitting HM version of this model with the same atmospheric component. This result is in contrast with the change from eddy-parameterized to eddy-permitting models, where there is a consistent increase in correlations and coupling Fig. 7 Wintertime climatology of ERA5 ∇ 2 SLP (top) ( 10 −9 Pa∕m 2 ) and NOAA OISSTv.2 ∇ 2 SST (bottom) ( 10 −10 • C∕m 2 ). The contours represent the NOAA OISST sea surface temperatures with contour interval 2 • C. The thick black line denotes the 18 • C isotherm coefficients among all models. However, this is consistent with other studies. Roberts et al. (2016), comparing an eddypermitting to an eddy-resolving version of the Met Office climate model, find no significant improvement of the SST and winds stress relationship. Bellucci et al. (2021), find little improvements on the relationship between the SST and THF when comparing an eddy-resolving configuration of the MPI-ESM model to the HighResMIP simulations of the same model.
Presently, we do not have a clear explanation for this behavior, but we note that the atmospheric resolution of 50 km is much coarser than the ocean resolution of 1/12 • ( ∼ 10km ) and this ocean eddy variability at those small scales cannot be seen by the atmosphere. If in the change from ocean eddy-permitting to eddy-resolving models a significant fraction of the ocean variability is transferred to those small scales, this will not be resolved by the atmosphere, resulting in weaker or similar correlations and coupling coefficients. Moreton et al. (2021), using three different configurations of the HadGEM3-GC3.1 model, they found evidence that increasing the oceanic resolution at a constant atmospheric resolution can deteriorate the representation of air-sea interaction. However, Bryan et al. (2010), comparing two eddy-resolving versions of the NCAR Community Climate System Model (CCSM3.5) with 0.5 • and 0.25 • atmospheric resolution, find no significant changes. They suggest that implementing more accurate parameterizations of the atmospheric boundary layer would be more beneficial than increasing the atmospheric resolution. To understand better these results, more in-depth analyses, in a multi-model approach, is required.
Despite the general tendency of the models to exhibit a better agreement with observations by enhancing the resolution, there is an inter-model spread which can be attributed to factors other than the spatial resolution. The coupling strength can be affected by the vertical resolution of the atmosphere and the number of vertical levels in the boundary layer (Chelton and Xie 2010). Moreover, different values for the winter period (DJF). The stippling shows statistically non-significant correlations at 95% significance level using two-sided student's t-test physical parameterizations in the convection and boundary layer schemes can significantly alter the surface wind response to the SST variability (Song et al. 2009;Perlin et al. 2014).
An important impact of the underestimation of the air-sea coupling that is not studied here is that on the storm track activity. Piazza et al. (2016) show that a simulation forced by smoothed SST fields over the GS region, leads to a northward shift of the storm track density towards the North America coast. Similar results are obtained by Woollings et al. (2010) and Small et al. (2014). This shift of the main storm track can also result in reduced storm activity over remote regions such as Greenland and the Mediterranean (Piazza et al. 2016). Haarsma et al. (2019) analysing coupled seasonal forecast experiments at different resolutions, find stronger injection of wave activity into the storm track with increasing oceanic resolution. In order to assess these impacts further, dedicated experiments have been conducted, and the results will be presented in future work.
Finally, even though we focus our analysis on the GS, the ocean forcing on the atmosphere has been also identified in other regions with strong ocean fronts and eddy activity, such as the Kuroshio and its Extention (KE) (e.g. Nonaka and Xie 2003;Putrasahan et al. 2013), the Antarctic Circumpolar Current (e.g. O'Neill et al. 2003), the Agulhas Return Current (e.g. O'Neill et al. 2005;Song et al. 2009) and the Brazil-Malvinas confluence (e.g. Tokinaga et al. 2005;Pezzi et al. 2009). Both observational and model analysis indicate that while the coupling coefficients are comparable for the KE and GS, higher values are retrieved in the Southern Ocean (Chelton et al. 2004;Maloney and Chelton 2006;Bryan et al. 2010). Despite the regional differences of the coupling magnitude, which might be related to differences in the MABL structure, we generally expect similar resolution dependence of the coupling strength in all the regions. Past studies have shown that in eddy-parameterized models the coupling is weak or even absent over all regions with intense mesoscale activity, indicating the oceanic resolution as a key feature for the air-sea coupling (e.g. Bryan et al. 2010;Small et al. 2019). However, we should be careful extrapolating our results to other regions since discrepancies may arise due to different intensities and orientations of the ocean fronts and/or air-sea temperature contrasts. For example, Bellucci et al. (2021) found that while eddy-parameterized models are able to reproduce the ocean forcing on the atmosphere over GS and KE, this is not the case for the Southern Ocean. This may be related to geographical dependencies of the transition length scales setting the transition from ocean-driven to atmosphere-driven regimes (Bishop et al. 2017;Small et al. 2019).

Conclusions
The GS region, with its sharp SST front and large eddy activity, is an area with strong ocean-atmosphere interactions, especially during winter when cold, dry continental air-masses are passing over. It is also a source region for atmospheric baroclinic instability, shaping the strength and position of the storm tracks that affect the weather and climate in Europe. Simulating accurately the ocean-atmosphere interaction is therefore key for correctly simulating the North Atlantic jet and storm track and European climate (Haarsma et al. 2019;Piazza et al. 2016;Small et al. 2014).
Using the simulations of six GCMs made according to the HighResMIP protocol, we have evaluated the ocean-atmosphere interaction and compared it with reanalyses and available observations. Two mechanisms have been investigated in detail: VMM and PAM. A clear dependence on resolution is found with increased ocean-atmosphere coupling strength and better agreement with reanalysis and observations, for both increasing ocean and atmosphere resolution, in agreement with Czaja et al. (2019). Based on the available simulations we conclude that eddy-permitting ocean resolution (0.4 • to 0.25 • ) and comparable atmosphere resolution are required for a realistic simulation of the two mechanisms. Configurations of eddy-parameterized ocean resolution exhibit lower correlation and coupling coefficients, indicating that the ocean-atmosphere coupling is weaker in the GS region. An intriguing result is that the only model with an eddy-resolving ocean revealed weaker or similar coupling and correlation coefficients compared to its eddy-permitting version. Further analyses is required to better understand this result.
Since most GCMs used to inform the IPCC AR5 and AR6 are eddy-parameterized models, we argue that their representation of the ocean-atmosphere interaction along the GS and other western boundary currents is flawed, with implications for the dynamics of the storm tracks and the associated weather and climate effects. Based on our results, we conjecture that they underestimate both the PAM and the VMM. Because of the limited length of the simulations and the many other drivers that affect the storm tracks when the resolution is changed (drivers that originate in the tropics or in the Arctic), we could not investigate how incorrectly representing ocean-atmosphere interaction affects the storm track dynamics. New, designed on purpose experiments are presently being analyzed.
acknowledge PRIMAVERA funding received from the European Commission under Grant Agreement 641727 of the Horizon 2020 research program.
Funding This work was supported by Horizon 2020 Framework Programme (641727).

Code availability
The computer code used in the present study is available from the corresponding author upon reasonable request.

Conflicts of interest
The authors declare that they have no conflict of interest.
Ethics approval Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.