North Atlantic gyre circulation in PRIMAVERA models

We study the impact of horizontal resolution in setting the North Atlantic gyre circulation and representing the ocean–atmosphere interactions that modulate the low-frequency variability in the region. Simulations from five state-of-the-art climate models performed at standard and high-resolution as part of the High-Resolution Model Inter-comparison Project (HighResMIP) were analysed. In some models, the resolution is enhanced in the atmospheric and oceanic components whereas, in some other models, the resolution is increased only in the atmosphere. Enhancing the horizontal resolution from non-eddy to eddy-permitting ocean produces stronger barotropic mass transports inside the subpolar and subtropical gyres. The first mode of inter-annual variability is associated with the North Atlantic Oscillation (NAO) in all the cases. The rapid ocean response to it consists of a shift in the position of the inter-gyre zone and it is better captured by the non-eddy models. The delayed ocean response consists of an intensification of the subpolar gyre (SPG) after around 3 years of a positive phase of NAO and it is better represented by the eddy-permitting oceans. A lagged relationship between the intensity of the SPG and the Atlantic Meridional Overturning Circulation (AMOC) is stronger in the cases of the non-eddy ocean. Then, the SPG is more tightly coupled to the AMOC in low-resolution models.


Introduction
The redistribution of heat in the North Atlantic Ocean is a key mechanism for the global climate system. Ocean circulation accounts for more than 30% of the global poleward heat transport accomplished by the ocean-atmosphere system (von der Haar and Oort 1973;Trenberth and Solomon 1994). In the Atlantic ocean, the meridional overturning circulation is the dominant contributor to the oceanic meridional heat transport for latitudes south of 50°N. Still, the horizontal circulation in the North Atlantic is efficient in redistributing heat over the subpolar latitudes (Tiedje et al. 2012;Dong and Sutton 2002). The North Atlantic barotropic circulation consists of a double gyre system: a cyclonic circulation over the subpolar basin (the subpolar gyre; SPG) and an anticyclonic circulation in the subtropics (subtropical gyre; STG). The SPG, approximately spanning the 45°N-65°N latitude range (Rhein et al. 2011), is a dominant large-scale feature of the surface circulation of the northwest Atlantic (Higginson et al. 2011). On the other hand, with warmer and saltier waters, the STG represents an enormous reservoir of heat and salt (Schmitz and McCartney 1993). As a branch of the STG, the Gulf Stream transports tropical and subtropical waters to the subpolar North Atlantic. Both the SPG and the STG contribute to the Atlantic Meridional Overturning Circulation (AMOC) and the SPG dynamics has also impacts on the Arctic sea-ice Renold et al. 2010). Hatun et al. (2005), for example, have studied the influence of the subpolar gyre dynamics on the thermohaline circulation and they conclude that salinity inflow to the Arctic basin is controlled by the dynamics of the SPG on inter-annual to inter-decadal time scales. Thus, the lowfrequency fluctuations in the North Atlantic horizontal gyre circulation can potentially affect future climate variability.
Being a region where the ocean-atmosphere interactions are intense, the North Atlantic Ocean variability is highly dominated by changes in the North Atlantic Oscillation (NAO; Deshayes and Frankignoul 2008;Robson et al. 2016), the main mode of inter-annual to inter-decadal 1 3 variability in the region. The SPG is influenced by the NAO through wind forcing and surface buoyancy fluxes (Eden and Willebrand 2001;Bellucci and Richards 2006;Böning et al. 2006;Deshayes and Frankignoul 2008). Observational and modelling studies suggest that the surface wind stress highly affects both, the strength and variability of the SPG (Curry et al. 1998;Böning et al. 2006;Häkkinen et al. 2011). The density structure is also an important component on the gyre transport (Mellor et al. 1982;Greatbatch et al. 1991;Eden and Willebrand 2001;Montoya et al. 2011). While buoyancy forcing is mainly responsible for the inter-decadal ocean variability, the wind stress changes act on shorter time scales (Eden and Jung 2001). From forced ocean simulations of the Atlantic, Eden and Willebrand (2001) found an instantaneous barotropic and a delayed baroclinic response of the ocean to a positive phase of NAO. The barotropic response consists of an anomalous anticyclonic circulation on the subpolar front. The baroclinic response leads to an intensification of the SPG. Accordingly, Marshall et al. (2001) have introduced the term 'inter-gyre gyre' as a circulation anomaly induced by Rossby wave adjustments to changes in the meridional drift of wind stress pattern produced by the NAO. In particular, the inter-gyre gyre would transport heat northward when the NAO is in its positive phase (Czaja and Marshall 2001;Bellucci et al. 2008).
The SPG cyclonic circulation also modulates the penetration of subtropical salty waters (Hatun et al. 2005) with implications in the variation of the mixed layer depth and consequently on water mass formation. A correlation between the SPG and the AMOC was also found in numerical models in association with increased deep convection (Delworth et al. 1993;Eden and Willebrand 2001;Spall 2008;Yoshimori et al. 2010;Born and Mignot 2011). At the same time, the SPG circulation could have two stable states, switching between a strong and a weak mode as a consequence of positive feedback mechanisms associated with the stratification (Born et al. 2013(Born et al. , 2015Born and Stocker 2014).
Observations on basin-scale are still limited when it comes to transports and process-related quantities. Estimating the mean circulation in the subpolar North Atlantic based on observations might be seasonally biased due to the sparse and quite short records of observations (Higginson et al. 2011). Therefore, ocean and climate models are useful tools in studying such processes. While the ability of global climate models (GCMs) in representing the AMOC has been regularly assessed across different generations of models (e.g. Cheng et al. 2013;Danabasoglu et al. 2014Danabasoglu et al. , 2016Roberts et al. 2020), relatively less attention has been devoted to the representation of the gyre circulation in stateof-the-art climate models. Since the dynamics of the North Atlantic have impacts on the distribution of Arctic sea-ice, the AMOC and the European climate, understanding how it is represented in the last-generation climate models is important to better simulate the low-frequency variability and therefore to improve the decadal climate predictions.
Thanks to the advances in computational power, considerable developments in climate models have been made in the last decades. A feature that was particularly recognized is the increased resolution. There is an amount of evidence that an enhanced resolution in climate models can improve the representation, not only of small-scale processes but also of the mean state of the climate. High-resolution eddy-permitting ocean models can adequately reproduce the salient features of the gyre circulation (Eden and Böning 2002;Treguier et al. 2005) and in particular, the Gulf Stream separation and the position of the North Atlantic Current (Marzocchi et al. 2015). However, the high-resolution in the ocean could be unimportant to improve the low-frequency variability of the gyre (Böning et al. 2006). It can also be expected that the enhanced resolution in the atmospheric component also could affect the oceanic variability because of the strong seaair interaction taking place in the North Atlantic.
In the context of the EU Horizon 2020 PRIMAVERA project, targeted climate simulations were run following the protocol of the High-Resolution Model Inter-comparison Project (HighResMIP; Haarsma et al. 2016). Several modelling groups have contributed with the same climate simulations in at least two configurations: a standard-resolution and a high-resolution. In some cases, the high-resolution set up consists of an increased resolution in both components, atmosphere and ocean. In other cases, the resolution was increased only in the atmosphere. In this paper, we take advantages of the PRIMAVERA models to investigate the sensitivity of the North Atlantic gyre circulation to the ocean and atmosphere horizontal resolution by using five different state-of-the-art coupled GCMs. We aim at assessing the added value of increased resolution in representing the mean state of the gyre circulation, the ocean-atmosphere interactions and the driven mechanisms dominating the lowfrequency variability in the North Atlantic Ocean.
In the next section, we will describe the model simulations and the methodology used for the analysis. The results regarding the mean state of the barotropic circulation and its variability will be presented in Sect. 3. We will discuss the results in Sect. 4 and summarize the conclusions in Sect. 5.

Models
In this study, we analyse historical simulations (hist-1950) covering the 1950-2014 period, performed following the Coupled Model Inter-comparison Project, Phase 6 (CMIP6) HighResMIP common protocol (Haarsma et al. 2016). The HighResMIP protocol has been specifically designed for the systematic investigation of the impact of horizontal resolution on the representation of climate-relevant processes in global circulation models. Five coupled models were selected, each one presenting a standard and a high-resolution configuration. For three of these models, EC-Earth3P (Haarsma et al. 2020), ECMWF-IFS (Roberts et al. 2018) and HadGEM3G-C31 , the resolution enhancement is achieved by refining both oceanic and atmospheric grids. In these models, the resolution of the ocean components ranges between non-eddying configuration (1° × 1°) and eddy-permitting (1/4° × 1/4°), whereas the atmospheric resolution ranges between 250 and 25 km. Two models feature an enhanced resolution only for the atmospheric component, ranging from 100 to 25 km. These are CMCC-CM2 (Cherchi et al. 2019) and MPI-ESM1-2 (Gutjahr et al. 2019). The oceanic resolution is 1/4° × 1/4° (0.4° × 0.4°) in CMCC-CM2 (MPI-ESM1-2). The ocean component of the climate models EC-Earth3P, ECMWF-IFS, HadGEM3-GC31 and CMCC-CM2 is based on the Nucleus for European Modelling of the Ocean framework (NEMO, Madec et al. 2016), whereas MPI-ESM1-2 uses the Max Planck Institute Ocean Model (MPIOM, Jungclaus et al. 2013). Some of the modelling groups have performed several ensemble members. When available, they are used to show the members' spread. For details on the models refer to Docquier et al. (2019), whereas the information regarding the simulations used in this study is provided in Table 1.
The model outputs analysed in this paper for EC-Earth3P

Methods
For the NEMO models, the North Atlantic barotropic streamfunction (BSF) and the meridional overturning streamfunction were computed with the diagnostic package CDFTOOLS (CDFtools 2020) using the 3D fields of horizontal velocity (u, v) and the meridional velocity (v), respectively. For MPI-ESM1-2, the model output of barotropic mass streamfunction and the meridional mass streamfunction were divided by the reference density of 1025 kg m −3 . The resulting BSF fields were remapped on a regular 1° × 1° grid to facilitate the inter-comparison. The meridional streamfunction in the Atlantic basin was used to compute the AMOC index. It was calculated as the maximum value of Atlantic streamfunction at 26.5°N and 53°N in the 900-1200 m depth interval.
The outputs of sea level pressure and wind stress from the atmospheric models were used to explore the sea-air interaction. The NAO index was computed as the leading Empirical Orthogonal Function of winter (DJFM) sea level pressure anomalies over the Atlantic sector, 20°N-80°N, 90°W-40°E (NCAR 2020). The wind stress curl was computed using the components of the wind stress (tauuo, tauvo) for each model.

The barotropic streamfunction
The mean BSF computed over the 1950-2014 period for high and low-resolution model configurations are shown in the first and second column of Fig. 1, respectively. To highlight the impact of the increased resolution, difference patterns between high and low resolution are also shown in the third column of Fig. 1. Dots indicate areas where the differences between the mean values (|HRmean − LRmean|) are higher than the square root of the two variances added up (sqrt(HRvariance + LRvariance)). We consider in those cases that the difference is statistically robust. Results from one single ensemble member for each model and configuration are shown. In all the cases, the climatology displays the subpolar and subtropical gyres (SPG and STG) as the minimum and maximum values of BSF, respectively ( Fig. 1). Largest and robust differences are found when both atmosphere and ocean resolutions are increased (upper three rows). In particular, the BSF fields are smoother with the lower resolution versions. Overall, the enhanced resolution is associated with a strengthening of the SPG cyclonic circulation (particularly marked along the western Greenland coastal boundary and in the Labrador basin), and a larger transport in the Gulf Stream (third column of Fig. 1). The position of the inter-gyre zone is around 45°N and it does not change with the resolution. We do not find systematic differences in the 2D gyre structure when the resolution is changed only in the atmosphere (last two rows of Fig. 1). Moreover, barotropic transports in the SPG and the Gulf Stream result slightly higher in the case of the lower atmospheric resolution for MPI-ESM1-2. We now focus on the mean strength of the northern cyclonic gyre by looking at the horizontal section outlined in the upper left panel of Fig. 1, which is similar to Results of five climate models in their high (first column) and low (second column) horizontal resolutions. The third column shows the difference between the two configurations. Dots highlight areas in which the differences are robust, that is, when the differences are larger than the square root of the variance of the two series summed up. One ensemble member for each model configuration is plotted the Overturning in the Subpolar North Atlantic Program (OSNAP) section. The climatological BSFs along the transect for each model configuration are plotted in Fig. 2. The dotted black line displays the barotropic streamfunction along the same transect as estimated from observations in Colin de Verdiere and Ollitrault (2016). In their paper, Colin de Verdiere and Ollitrault (2016) estimate the barotropic circulation of the global ocean under geostrophic and mass conservation assumptions using Argo float displacements and the World Ocean Atlas (WOA) climatology. For the estimates shown in Fig. 2, Argo floats data till 2015 were also included (Ollitrault et al., 2020). Models with different horizontal resolutions in the atmosphere and the ocean are plotted in Fig. 2a whereas models with different resolution only in the atmosphere are plotted in Fig. 2b. Solid (dashed) lines display the high (low) resolution versions. When more than one ensemble member is available, the line represents the ensemble mean and the dashed area illustrates the members' spread. The spread among the ensemble members is low indicating that the climatology is similarly represented by each member. In general, the maximum barotropic transports in the SPG (absolute value of BSF) are found along the western section, around 56°W and over the eastern coast of Greenland, whereas smaller absolute values of BSF are found in the eastern part of the transect. It is found that the impact of resolution does generally overcome the model-tomodel differences when both atmosphere and ocean resolutions are increased (Fig. 2a). In fact, the high-resolution versions of EC-Earth3P, ECMWF-IFS, HadGEM3-GC31 perform very similarly. The lower resolution versions simulate smaller transports all along the transect. It is worth noting, however, that the ocean model component is NEMO in all those cases, although in somewhat different configurations. The transports simulated by the non-eddying configurations are closer to the observational estimates than the ones simulated by the eddy-permitting models. However, the WOA temperature and salinity climatology is used to estimate the observational barotropic transports and therefore they might be underestimated. On the other hand, when the resolution is increased only in the atmosphere, differences  Fig. 1 for a models with different horizontal resolutions in both components, and b models with different horizontal resolution in the atmosphere. The lines represent the ensemble mean and the shaded area the members' spread. The number of ensemble members analysed for each model configuration is indicated within brackets. The dotted black line indicates the barotropic transport estimated from observations among models are larger than the differences implied by the enhanced resolution (Fig. 2b). In particular, compared to CMCC-CM2, MPI-ESM1-2 simulates lower barotropic transports in the western tongue and less spatial variability in the eastern one.

SPG and STG indices
As stated in the previous section, the maximum absolute value of BSF across the SPG (Fig. 2) is reached, in general, near the western end of the section. Following Danabasoglu et al. (2014), we computed the SPG index as the minimum BSF between 65° and 40°W at 53°N; while the STG index is derived as the maximum BSF between 80° and 60°W at 34°N. The time series of the annual SPG and STG indices are plotted in Fig. 3. The SPG index ranges roughly between 20 and 35 (35-55) Sverdrup (Sv; 1 km 3 s −1 ) for the low (high) version when both, atmosphere and ocean resolutions are increased (Fig. 3a). For CMCC-CM2 and MPI-ESM1-2 values range between 45 and 75 (30-50) Sv, respectively (Fig. 3b). Xu et al. (2013), based on observational data from Fischer et al. (2004) and Fischer et al. (2010), estimated the southward transport in the Labrador Sea at 53°N to be 37-42 Sv. Therefore, the high-resolution models shown in Fig. 3a and the transport simulated by MPI-ESM1-2 are within the estimated range. The low-resolution versions in Fig. 3a underestimate it whereas CMCC-CM2 in its two configurations is the only model that overestimates the observed transports. No long-term trend in the SPG index is evident (Fig. 3a, b), except for CMCC-CM2-HR4 in which a decrease of almost 30 Sv is simulated within the 1986-2007 period.
The STG index ranges between 25-40 (55-95) Sv in the cases of low (high) resolution models plotted in Fig. 3c, and between 65-105 (50-75) Sv for MPI-ESM1-2 (CMCC-CM2; Fig. 3d). All the eddy-permitting oceans present higher inter-annual variability in the STG index compared to their low-resolution counterparts, probably explained by a more realistic representation of mesoscale eddies associated with the Gulf Stream. Curiously, MPI-ESM1-2 is the only model that simulates higher values of the SPG and STG indices in its low-resolution version. No long-term trend is evident in the STG index within the entire simulated period, from 1950 to 2014.
In general, there is no clear relationship between the SPG and STG indices for the same model setup (Fig. 4a, b), except for MPI-ESM1-2-HR, in which a slightly positive relationship is found. However, models with coarse ocean and atmosphere resolutions reproduce the minimum transports in both gyres (Fig. 4a). Again, the differences between the models that change the resolution only in the atmosphere (CMCC-CM2 and MPI-ESM1-2) are higher than the differences due to the increased resolution (Fig. 4b). The relationship between the annual SPG index and the annual AMOC index at 26.5°N is plotted in Fig. 4c, d. High values of SPG index are always associated with high values of the AMOC and vice versa. That is, models that simulate higher barotropic transports inside the SPG, are also able to reproduce a stronger AMOC. This is particularly evident when ocean resolution is enhanced.

The 2D-fields of BSF
We find no systematic long-term trend in the BSF evolution, either among the various models or amid the two configurations of a single model (Fig. 3). Only the CMCC-CM2 model in its low-resolution displays a relatively strong weakening of the SPG (Fig. 3b).
To investigate the inter-annual variability of the BSF, we performed an Empirical Orthogonal Function (EOF) analysis in the North Atlantic domain within 260-360ºE and 20-70ºN. We restrict our analysis to the 4 leading EOFs. The first EOF explains around 30% of the variance in the non-eddy ocean models (EC-Earth3P, ECMWF-IFS-LR and HadGEM3-GC31-LL) and in the two configurations of MPI-ESM1-2 (0.4º × 0.4º ocean resolution); and between 10 and 15% in the eddy-permitting ocean models (EC-Earth3P-HR, ECMWF-IFS-HR, HadGEM3-GC31-HM and CMCC-CM2). The EOF1 patterns are shown in Fig. 5, where only one ensemble member for each model is displayed. The first mode of variability generally consists of a south-west to north-east oriented anomaly of the circulation in the middle of the North Atlantic, centred at ~ 40ºN or slightly south.
The pattern is similar in all the cases, although less marked in CMCC-CM2 models, and it is consistent with variations in the position of the inter-gyre zone Bellucci et al. 2008). This is more evident in the cases of non-eddy oceans, in which the explained variance associated with this mode is even higher since the eddy-permitting models allow for smaller-scale variability. The time series of the corresponding principal component (PC1) and the relative power spectral density (PSD) are shown in Fig. 6. For the PSDs, the thick curves are the spectra of the PC1 time series whereas the thin lines represent the 90% confidence level based on the theoretical Markov spectrum for each series. To take into consideration all the ensemble members (when more than one is available), the detrended PC1 time series were concatenated to obtain one single spectrum for each model configuration (Fig. 6e). The low-frequency variability associated with EOF1 peaks around 11 years for the high-resolution configurations of ECMWF-IFS, HadGEM3-GC31 and EC-Earth3P (although not significant, Fig. 6c) and the low-resolution version of MPI-ESM1-2 (Fig. 6d). A significant peak around 6 years is found for EC-Earth3P-HR (Fig. 6c) and around 5 years for CMCC-CM2-HR4 (Fig. 6d). EOF1 patterns and the corresponding explained variances are very similar among the ensemble members (not shown). Conversely, the PSDs when all the ensemble members are considered show a relatively high spread. However, the peak around 11 years is still present for the high-resolution configurations of ECWMW-IFS, HadGEM-GC31 and EC-Earth3P (although this is only significant for EC-Earth3P) and the peak around 6 years in EC-Earth3P-HR is significant even if all the ensemble members are considered.  In summary, the first mode of inter-annual variability in the BSF is related to variations in the position of the intergyre zone. The variance explained by this mode is higher in the cases of the non-eddy models compared with the eddypermitting models. However, the energy associated with the 5-11 years variability of the PC1 is generally higher in the high resolutions in the ocean. To explore the forcing mechanisms of the first mode of variability, the wind stress curl for each model is correlated to the PC1s time-series (Fig. 7). In most cases, the highest negative correlations are located at ~ 45°N or a few degrees south. More generally, the patterns of correlation between PC1 and the wind stress curl are similar to those of the EOF1 for the BSF (Fig. 5): a negative (positive) anomaly in the wind stress curl (Fig. 7) is interpreted as a signature of the wind stress forcing an anomalous anticyclonic (cyclonic) barotropic circulation. Therefore, the first mode of inter-annual variability (5-11 years) of the BSF is generally driven by wind and represents variations in the inter-gyre position. The relation between the wind-stress and the barotropic transports is particularly strong in both configurations of EC-Earth3P, ECMWF-IFS and MPI-ESM1-2 and HadGEM3GC31-LL.

Relation with the North Atlantic Oscillation (NAO)
The correlation patterns between the NAO index and the wind stress curl are shown in Fig. 8. In all models and configurations, the wind stress response to the NAO is similar, with positive anomalies of wind curl along the Greenland coast, negative anomalies at midlatitudes (30°-50°N) and again positive anomalies south of 30°N (Fig. 8). This pattern is obtained independently from the atmospheric resolution. Except for HadGEM3-GC31-HM and CMCC-CM2, this correlation pattern generally resembles the positive phase of the correlation pattern between wind curl stress and PC1 (Fig. 7). Specifically, on these times scales, the wind stress curl pattern associated with NAO + leads to a coincident northward shift in the inter-gyre gyre. This would appear to be consistent with a negative (positive) anomaly of the wind stress curl (barotropic streamfunction) around 45°N. Accordingly, the combination of the results in Figs. 7 and 8 suggests that the first mode of inter-annual variability of the BSF in the North Atlantic is typically forced by the NAO through an ocean barotropic response to the wind stress.
To assess this forcing mechanism in the different models and to explain the different behaviour in HadGEM3-GC31-HM and CMCC-CM2, we analysed the cross-correlations between PC1 and NAO index (Fig. 9). In all simulations, positive and significant correlations are obtained for time lags between 0 and 5 years, NAO leading PC1. More specifically, the correlation generally peaks at about 0 lag with the sole exceptions of HadGEM3-GC31-HM and CMCC-CM2, for which the maximum correlation is obtained at time Fig. 7 Patterns of correlation between the first principal components PC1 and the wind stress curl for each model configuration. One ensemble member for each model configuration is plotted lags between 3 and 5 years, NAO leading PC1. This could explain the relatively low correlations obtained between PC1 and the wind stress curl (Fig. 7), suggesting that in these cases a delayed baroclinic response of the gyre circulation to the NAO is involved (Anderson et al. 1979;Bellucci and Richards 2006).
The previous analysis showed that, in the majority of the inspected models, the barotropic circulation in the North Atlantic instantaneously responds to the NAO by changes in the inter-gyre position. We computed the cross-correlation between the SPG index and the NAO index to identify a response in the intensity of the gyre. The time-series were filtered to retain the low-frequency variability before computing the cross-correlations, and a 7 years cutoff period was used for the low-pass filter. The resulting lagged correlations are shown in Fig. 10, where the NAO index leads SPG index. One ensemble member for each model configuration is plotted in Fig. 10a and the mean correlation (computed as the average across individual members' correlations) together with the ensemble' spread are plotted in Fig. 10c. For all numerical systems that include the NEMO model at highresolution, the correlation between SPG and NAO indices is maximum, positive and statistically significant when NAO leads SPG for about 2-4 years. Similar results are obtained in the other cases except for MPI-ESM1-2-XR, although correlations are not significant. This result implies that a relatively maximum barotropic transport inside the SPG is obtained around 3 years after a positive phase of NAO. Therefore, the response of the barotropic circulation in the North Atlantic to the NAO seems to be: (a) an instantaneous variation in the position of the inter-gyre and; (b) a delayed intensification of the SPG. This is a common response in most of the models. While the instantaneous response is more evident in the non-eddy models, the delayed response is stronger in the eddy-permitting ones.

Relation with the Atlantic meridional overturning circulation (AMOC)
The overturning circulation relates to the formation of dense water in the centre of the subpolar gyre and within the Greenland, Iceland and Norwegian Seas, but it is still unclear which is the main driver setting the AMOC amplitude. The cross-correlations between the SPG index and the AMOC index at 53°N (subpolar gyre) and 26.5°N are presented in Fig. 11. Results of this analysis are model and member dependent and the correlations are not always significant. One ensemble member for each model configuration is plotted in Fig. 11. For the significance test, the 95% confidence level was calculated as 2/sqrt(N), where N is the number of independent data based on the e-folding time-scale of the autocorrelation. In general, low-resolution configurations present the maximum positive correlations (although not always significant) with the AMOC at

Discussion
We investigate the impacts of the enhanced horizontal resolution in the ocean and atmosphere on the North Atlantic barotropic circulation by analysing numerical simulations from five different state-of-the-art global climate models following the HighResMIP protocol. Independently from the modelling system and the horizontal mesh, the structure of BSF climatologies is consistently Fig. 9 Cross-correlation between the NAO index (leading) and the PC1. Series were filtered with a 7 years low-pass filter before computing the cross-correlations. The area enclosed by dotted lines rep-resents the 95% confidence levels calculated as 2/sqrt(N), where N is the number of independent data based on the time that takes autocorrelation to fall below 1/e. Positive lags mean NAO leads PC1 well reproduced in comparison with what is reported in the literature. The inter-gyre position is generally around 45°N in all the considered runs as found, for instance, by Born et al. (2013) who reported the separation of the SPG and the STG between 45°N and 50°N in the CMIP3 simulations. We obtain a strength of the SPG ranging from 20 to 55 Sv. These values are consistent with the range 25-60 Sv reported by Treguier et al. (2005) in the western edge of the SPG from an ensemble of regional models. Differences in the structure of the BSF are found when both ocean and atmosphere resolutions increase but no systematic differences are found if only the atmospheric resolution changes. Therefore, we might infer that enhanced resolution in the ocean is responsible for the changes in the climatological field of BSF, particularly when this is increased from non-eddy to eddy-permitting. Apart from the small scale features due to the presence of mesoscale eddies, the eddy-permitting oceans display the maximum SPG transport to the western edge of the gyre. This is consistent with previous studies reporting the SPG centred to the east in coarser ocean models (Born et al. 2013). In general, the eddy-permitting ocean models exhibit stronger transports inside both gyres and are also able to simulate a stronger meridional overturning in agreement with other studies (Roberts et al. 2020). For some of the model couplets (notably, those changing the resolution of both the atmospheric and oceanic components) several ensemble members were analysed, for each (low and high-resolution) model configuration. Regarding the climatological values and the inter-annual variability of the SPG index, the ensembles' spread is in all those cases much smaller than the intra-model differences, indicating that our main results are not sensitive to the ensemble size. The response of the climate system to resolution is different when only the atmosphere grid is refined. In fact, in the MPI-ESM1-2 simulations, the AMOC circulation is less intense with finer atmospheric resolution. A higher resolution in the atmospheric component of coupled models can better Fig. 10 Cross-correlation between the NAO index (leading) and the SPG index after filtering the time series with a 7 years low-pass filter. One ensemble member for each model configuration is plotted in a and the mean curve and members' spread are plotted in c. The area enclosed by dotted lines represents the 95% confidence levels calculated as 2/sqrt(N), where N is the number of independent data based on the time that takes autocorrelation to fall below 1/e. Positive lags mean NAO leads SPG, negative lags mean SPG leads NAO reproduce cyclones, yielding more chaos and therefore reducing the intensity of the mean winds (Sein et al. 2018). As a consequence, Sein et al. (2018) have found weaker AMOC in response to increased atmospheric horizontal resolution. The opposite occurs in stand-alone ocean simulations (Jung et al. 2014), which are forced with reanalysis that assimilates data independently of their resolution. We have found that the first mode of inter-annual variability explains more variance in the models featuring a lowresolution ocean grid. Both, the first EOF pattern and its explained variance feature a small intra-members' spread, so our conclusions regarding the comparison between different ocean model grids are weakly sensitive to the ensemble size. However, when looking at the dominant spectral peaks in PSDs of the first principal components, the spread among the ensemble members is relatively high. By comparing simulations with an eddy-resolving ocean (1/12° × 1/12°) against a lower-resolution one (1/3° × 1/3°), Böning et al. (2006) concluded that a higher resolution model is needed to realistically represent the mesoscale eddies, but it is not relevant for capturing the low-frequency variability of the gyre.
The relationship between the leading mode of inter-annual gyre variability with the NAO has also been inspected. In general, the barotropic circulation responds to the NAO in two ways. On the one hand, there is a rapid response to wind stress by shifting the position of the inter-gyre zone. This result is consistent with Eden and Willebrand (2001) who stated that changes in the wind stress produced by a positive NAO are associated with an anomalous anticyclonic circulation in the subpolar front at lag 0. On the other hand, there is an oceanic slow response by increasing the intensity of the SPG. Deshayes and Frankignoul (2008) performed a 50-year hind-cast simulation in the North Atlantic with a regional ocean set up to study the inter-annual to interdecadal variability in the oceanic circulation. In agreement with our results, the authors found that the SPG responds with a time lag of 3 years to a positive phase of NAO. Likewise, the SPG has been shown to undergo an intensification in 2-3 years following a positive phase of NAO, in several model studies (Häkkinen 1999;Eden and Willebrand 2001;Gulev et al. 2003).
Analysing the relationship between the SPG and the AMOC, we have found that the low-resolution models work better in capturing a lagged relation, SPG leading AMOC. This is also consistent with previous works. Yeager (2015) suggested that the relation between AMOC and SPG is Fig. 11 Cross-correlation between the AMOC index (leading) at a, b 53ºN and c, d 26.5ºN and the SPG index after filtering the time series with a 7 years low-pass filter. One ensemble member for each model configuration is shown. The area enclosed by dotted lines represents the 95% confidence levels calculated as 2/sqrt(N), where N is the number of independent data based on the time that takes autocorrelation to fall below 1/e. Positive lags mean AMOC leads SPG, negative lags mean SPG leads AMOC associated with the bottom pressure torque being responsible for decadal, buoyancy-forced changes in the gyre circulation, providing coupling between AMOC and SPG. Performing a multi-regional-model analysis, Danabasoglu et al. (2016) suggested that the SPG sea surface height (SSH) changes tend to lead AMOC changes by 2-3 years. This result implies that variations in the SPG could be used as a proxy for changes in the AMOC in the past , or as potential indicator of future changes in the AMOC (Yong-Qi and Yu 2008).
It is worth highlighting that four over five high-resolution ocean configurations analysed in this study, share the same ocean component, NEMOv3.6 code. Although either the sea-ice model or physical schemes and parametrizations (e.g. eddy parametrizations, horizontal momentum diffusion, background vertical eddy viscosity, background vertical eddy diffusivity, isopycnal tracer diffusivity, eddy-induced velocity coefficient) differ, this might constitute a limitation of this inter-comparison analysis.

Conclusions
By analysing historical simulations in the period 1950-2014 in five state-of-the-art climate models in two different horizontal resolutions, we yield the following conclusions: • The mean gyre circulation exhibits the largest differences when the resolution in the ocean is increased from noneddy to eddy-permitting. A stronger mean circulation of the double gyre system is found when both ocean and atmosphere resolutions are increased. However, no systematic changes were found when only the atmospheric grid resolution is modified. • There is some evidence that increasing resolution in the atmosphere alone determines a weaker AMOC, at least around 26.5°N. • Most of the examined models show no long-term trend in the gyre circulation strength during the analysed period. • The first mode of inter-annual variability appears to be primarily driven by the NAO in all the examined simulations. Non-eddy ocean models show an instantaneous barotropic response to the NAO through the wind stress forcing. This response manifests itself through the emergence of an inter-gyre gyre circulation, straddling the cross-gyre climatological boundary. Conversely, models including an eddy-permitting ocean feature a slow response to the NAO driven by buoyancy forcing, followed by an intensification of the SPG circulation. • The relationship between the SPG and the AMOC indices is weak in the analysed cases. However, the non-eddy models exhibit similar behaviour and are consistent with what is reported in the literature. In this sense, the SPG index is more tightly coupled to the AMOC in the lowerresolution models.