Understanding the sensitivity of the North Atlantic subpolar overturning in different resolution versions of HadGEM3-GC3.1

The Atlantic Meridional Overturning Circulation (AMOC) is a key component of the global climate but is not simulated consistently across models or model resolutions. Here, we use a hierarchy of the global coupled model HadGEM3-GC3.1, with ocean resolutions of 1°, ¼°, and 1/12°, to evaluate the subpolar AMOC and its sensitivity to horizontal resolution. In line with observations, the models show that the mean overturning and surface forced water mass transformation (SFWMT) are concentrated in the eastern subpolar gyre rather than in the Labrador Sea. However, the magnitude of the overturning along the OSNAP line at medium and high resolutions is 25% and 40% larger than in the observations, respectively. This disagreement in overturning strength is noted for both OSNAP


Introduction
The Atlantic Meridional Overturning Circulation (AMOC) is a key component of the global climate through its role in the transport of heat and carbon northward in the Atlantic (Lozier et al., 2010).Changes in the strength of the AMOC are thought to be related to a wide range of local and global climate impacts (Bellomo et al., 2021;Zhang et al., 2019).The AMOC has been studied over the past decades to understand its impact on the wider climate and to predict future changes (Jackson et al., 2022).In particular, coupled climate models are used to investigate the forced response of the AMOC to historical forcing over a large range of timescales.In a context of global warming, climate models predict a risk of substantial AMOC weakening in the coming years partly due to a shutdown of open-ocean deep convection in the Labrador and Greenland seas (Brodeau & Koenigk, 2016;Heuzé, 2017).However, recent observations showed that the AMOC weakening might be overestimated in climate models because the deep-convection in the Labrador Sea is not the primary observed source of dense water formed in the subpolar North Atlantic (Fu et al., 2020;Roquet & Wunsch, 2022;Thomas et al., 2015;Zhang & Thomas, 2021).On the contrary, the dense water is mainly formed over the eastern subpolar gyre (Chafik & Rossby, 2019;Lozier Abstract The Atlantic Meridional Overturning Circulation (AMOC) is a key component of the global climate but is not simulated consistently across models or model resolutions.Here, we use a hierarchy of the global coupled model HadGEM3-GC3.1,with ocean resolutions of 1°, ¼°, and 1/12°, to evaluate the subpolar AMOC and its sensitivity to horizontal resolution.In line with observations, the models show that the mean overturning and surface forced water mass transformation (SFWMT) are concentrated in the eastern subpolar gyre rather than in the Labrador Sea.However, the magnitude of the overturning along the OSNAP line at medium and high resolutions is 25% and 40% larger than in the observations, respectively.This disagreement in overturning strength is noted for both OSNAP East and OSNAP West, and is mainly due to anomalously large SFWMT rather than anomalously large interior mixing or overflow transport from the Nordic Seas.Over the Labrador Sea, the intensification of SFWMT with resolution is explained by a combination of two main biases.Anomalously warm surface water enhances heat loss and reduces the extension of marginal sea ice, which increases the surface density flux over the boundary of the basin.A bias in salinity leads to anomalously dense surface water that shifts the outcropping area of the AMOC isopycnal and results in intense dense water formation along the boundary of the basin at medium and high resolutions.Thus, our analysis sheds light on a range of model biases responsible for large overturning over the Labrador Sea in climate models.

Plain Language Summary
The circulation and formation of dense water over the North Atlantic plays a key role for the Northern Hemisphere climate as it influences sea surface temperature, the formation of hurricane and the intensity of rainfall over Europe.It is thus important to evaluate how this circulation will change in a context of global climate change.However, the simulated overturning is not necessarily consistent with observations and differs across models or model resolutions.In this work, we evaluate the impact of horizontal resolution for the overturning over the subpolar gyre in the climate model HadGEM3-GC3.1.We find that its magnitude at medium and high resolutions is 25% and 40% larger than in the observations, respectively.The water mass transformation induced by the air-sea exchanges (SFWMT) over the subpolar gyre mainly explains this disagreement.Over the Labrador Sea, the relatively strong SFWMT in these models is related to anomalously small sea ice coverage as well as anomalously warm and salty surface water.
PETIT ET AL. 10.1029/2023JC019672 2 of 16 et al., 2019) and the buoyancy forcing is identified as a main driver of the dense water formation over the subpolar gyre (Petit et al., 2020).This is reinforced by Jackson and Petit (2022) who recently showed that the surface forced water mass transformation (SFWMT) is a good indicator of the mean overturning over the subpolar gyre in CMIP6 models.Then, given our new understanding of the mechanisms driving the dense water formation in the subpolar North Atlantic, how can we explain the larger overturning over the Labrador Sea in climate models as compared to observations?
The representation of a realistic AMOC, which influences its forced response (Bellomo et al., 2021;Jackson et al., 2020;Roberts et al., 2020), remains a challenge in climate models.The representation of the AMOC can be affected by biased processes; for instance, the overflow transport from the Nordic Seas (Yeager & Danabasoglu, 2012), the exchanges between the boundary and interior of the basins (Georgiou et al., 2019(Georgiou et al., , 2020;;Spall, 2011;Tagklis et al., 2020), the deep convection (Danabasoglu et al., 2014;Koenigk et al., 2021;Oldenburg et al., 2021), the representation of mesoscale eddies (Moreton et al., 2020(Moreton et al., , 2021)), or the ocean circulation between the subtropical and subpolar gyres (Jackson et al., 2020).The misrepresentation of narrow boundary currents in low-resolution models can also affect the AMOC through the advection of anomalously warm and salty mode water (Talandier et al., 2014;Treguier et al., 2005).
Although these biases are expected to be overall stronger in models at low-resolution, it remains unclear if the representation of the AMOC is improved in higher-resolution models.In the CESM1 model (Chang et al., 2020;Zhang et al., 2020), the simulated overturning and SFWMT over the subpolar gyre are closer to observations in models at high resolution as compared to the low resolution (Oldenburg et al., 2022;Yeager et al., 2021).In particular, they found a strong overturning over the Labrador Sea at low resolution, while that at higher resolution was weaker and closer to the observations.However, a comparison of the AMOC at 26.5°N between a hierarchy of the HadGEM3-GC3.1 model at different resolutions revealed that the strong overturning in the experiments at high resolution are not closer to observations than the weak overturning in the experiments at low resolution (Roberts et al., 2020).Because there has not been a systematic study of the AMOC sensitivity over the subpolar gyre in HadGEM3-GC3.1, it is not clear if the differences between these studies are related to the models or to a meridional incoherence of the AMOC sensitivity between the subpolar and subtropical gyres.Indeed, the subpolar overturning estimated in HadGEM3-GC3.1 at a 1° and ¼° resolutions by Jackson and Petit (2022) showed a link between the salinity biases over the Labrador Sea and the AMOC, but the impact of resolution for the representation of the AMOC remains unclear (Menary et al., 2020).
Based on these papers, we investigate the resolution dependence of the subpolar overturning in HadGEM3-GC3.1 by extending the hierarchy of resolution used in our study, ranging from 1° to 1/12° for the ocean grid and from 250 to 50 km for the atmosphere grid.Because the OSNAP array provides an exceptional benchmark for the quantification of the AMOC over the subpolar gyre, we use the 47-months timeseries of the overturning at OSNAP to evaluate the simulated overturning in these models.The models and observational data sets are detailed in Section 2. The sensitivity of the dense water formation at OSNAP to model resolution is analyzed in Section 3, and processes explaining the differences are identified in Section 4. In Section 5, we investigate the role of air-sea fluxes and sea ice melting for the representation of dense water formation over the Labrador Sea.Sections 6 and 7 summarize and discuss the wider implications of these results.

Model HadGEM3-GC3.1
We use the 1950s control run of 3 simulations of the coupled model HadGEM3-GC3.1 (Roberts et al., 2019) carried out with the High Resolution Model Intercomparison Project (HighResMIP) protocol for CMIP6 (Haarsma et al., 2016).HadGEM3-GC3.1 is based on the Unified Model for the atmosphere, CICE for the sea ice and NEMO for the ocean (Roberts et al., 2020, and references therein).From low to high resolutions, we use HadGEM3-GC3.1-LLwith an ocean resolution of 100 km and an atmosphere resolution of 250 km, HadGEM3-GC3.1-MMwith an ocean resolution of 25 km and an atmosphere resolution of 100 km, and HadGEM3-GC3.1-HHwith an ocean resolution of 8 km and an atmosphere resolution of 50 km.We will refer to these as the LL, MM and HH models, respectively.Because it is difficult to cleanly separate the impacts of changes in the atmosphere and ocean resolutions, here we consider their combined effects on oceanic biases.While we do not further discuss these complex coupled interactions in the study, it is important to note that 10.1029/2023JC019672 3 of 16 oceanic biases can be caused by atmospheric effects (e.g., shift in jet stream due to resolution) or coupled effects (e.g., air-sea feedbacks, see Moreton et al., 2021).
The three models have the same vertical resolutions, with 75 levels in the ocean and similar parametrizations.The few parameters adjusted because of the change in horizontal resolution are described in Section 2.2 of Roberts et al. (2019).The higher resolutions MM and HH are typically described as eddy-permitting and eddy-rich resolutions because they partly resolve mesoscale eddies up to the mid-latitudes (∼30°N/30°S) and high-latitudes (∼60°N/60°S), respectively (Hallberg, 2013;Moreton et al., 2020).Hence, they do not include a parametrization of eddy transport.On the contrary, the lower resolution LL uses standard eddy parametrizations (Gent & Mcwilliams, 1990;Kuhlbrodt et al., 2018;Redi, 1982).Because of the dependence of the wind stress on surface currents, it is also expected that its representation will differ between eddy-parametrized and eddy-rich models (Hewitt et al., 2020).
The model version differs slightly from those used by Menary et al. (2020) because the HighResMIP simulations use the simplified aerosol forcing from the second version of the Max Planck Institute Aerosol Climatology (MACv2-SP) scheme (Stevens et al., 2017), rather than the Global Model of Aerosol Processes (GLOMAP-mode) scheme (Mulcahy et al., 2018).Boundary conditions (e.g., atmospheric composition) are fixed and use constant 1950s forcing.Note that comparing simulations using stabilized external forcing (i.e., control simulations) with observations can lead to uncertainties (King et al., 2020), although they cannot be robustly quantified in our study.Because the HH simulation is only 100-years long, we use only the first 100 years from each simulation for consistency.A short spinup of 30 years initializes the simulations, which is not included in the 100 years of simulation analyzed here.More details on the simulations can be found in Roberts et al. (2019).

Observations
The OSNAP mooring array provides continuous measurements of volume, heat and freshwater transports through the subpolar North-Atlantic.It is composed of the OSNAP East section from Scotland to the southeast tip of Greenland, and the OSNAP West section from the southwest tip of Greenland to the Labrador shelf (Figure 1a).These 30-days mean fields are gridded with a horizontal resolution of ∼25 km and a vertical resolution of ∼20 m (Lozier et al., 2019).The overturning at OSNAP is estimated from the velocity fields across each section (Li et al., 2017) during 3 full years, from January 2015 to December 2017.Note that this time series is longer than the 21-months timeseries used by Menary et al. (2020).
As described by Petit et al. (2020), we consider the southward transport of overflow water across the GSR of 6.6 ± 0.4 Sv.This value has been estimated by combining the southward overflow transport at the Faroe-Bank Channel (Hansen et al., 2016) with that at the Denmark Strait (Jochumsen et al., 2017).These transports are provided as monthly means by the Atlantos project from 1996 to 2017 (de Young et al., 2019).An additional overflow transport of 1.2 Sv was included to these overflow transports to account for Iceland-Scotland Overflow Water inflow through the Iceland-Faroe Ridge (Sherwin, 2010) and the Wyville-Thomson Ridge (Hansen et al., 2018).Further details on the specific acquisition of each of these observations can be found in the referenced publications.

The Ocean Reanalysis ECCO
The simulated surface forced water mass transformations (SFWMT) are evaluated in comparison to the SFWMT estimated from ECCOv4r4 (ECCO).The monthly heat and freshwater fluxes and the sea surface salinity and temperature fields are gridded onto a 0.5° grid from 1992 to 2017 (Forget et al., 2015).We consider ECCO as an analog to observations for SFWMT because the influence of sea ice is included in its freshwater fluxes.In contrast, computations of SFWMT from atmospheric reanalysis (e.g., NCEP or ERA5) only include evaporation and precipitation in their estimations of freshwater fluxes.Moreover, previous studies documented the consistency between the observed overturning at RAPID and OSNAP with that estimated from ECCO (Jackson et al., 2019).

Estimation of the Overturning and Surface Forced Water Mass Transformation
The simulated overturning is estimated in density space along the two OSNAP sections and along the Greenland-Scotland Ridge (GSR) for the 3 models (Figure 1).The sections are constructed to extract velocities as close as possible to the observation lines while conserving the total model transport (Figure S1 in Supporting Information S1).The overturning stream function is estimated by summing the transports along the section and by accumulating them from the bottom to the surface.The overturning stream function is calculated in density space with density referenced to the surface, σ 0 , for comparison with the OSNAP observations.The stream function is computed using monthly mean data and averaged over the 100 years of the model experiments to form the model climatology.All uncertainties are expressed as the standard deviation of their interannual variability.As detailed in Menary et al. (2020), no compensation due to the net throughflow is applied in our estimation of overturning at OSNAP.We instead integrate transports from the bottom upwards to better compare the deep layers of the models with observations.Following Menary et al. (2020), the overturning across each section is defined as the maximum of the averaged overturning stream functions shown in Figure 1, which avoid potential artifacts due to changing the density of the stream function maximum each month.The density associated with this maximum is called the MOC isopycnal,   moc .The overturning divergence over the area bounded by the OSNAP East and GSR sections, called Irminger and Iceland Basins, is defined as the maximum of the difference between the overturning stream functions at these two sections.
The dense water formation is partly forced by surface fluxes, which can be quantified with the SFWMT assuming a steady state.As described in previous works (Desbruyères et al., 2019;Josey et al., 2009;Petit et al., 2020;Speer & Tziperman, 1992;Walin, 1982), the SFWMT is estimated across each isopycnal σ by integrating the buoyancy forcing (B) over the surface area of each density bin, Δσ: 10.1029/2023JC019672 5 of 16 where, In Equation 1, C p is the specific heat, α is the thermal expansion coefficient and β is the haline contraction coefficient.The surface heat flux (Q) includes radiative, turbulent, latent, frazil and heat content fluxes as the "hfds" output of the models.The surface freshwater flux (   net ) is estimated by combining precipitation ("pr" and "prsn"), evaporation ("evs"), river runoff ("friver") and ice processes ("ficeberg" and "fsitherm").The surface density is estimated for each month from the temperature (T) and salinity (S) fields at surface, such that the SFWMT is set to zero when σ does not outcrop in a given month over the considered area.Here, we use a density bin of 0.1 kg m −3 and we define two areas (Figure 1a) over which the SFWMT is integrated: over the Labrador Sea (LS) between 67°N and OSNAP West (SFWMT LS ) and over the Irminger and Iceland Basins (IIB) between OSNAP East and GSR (SFWMT IIB ).Following Petit et al. (2020), we consider that the dense water formation (i.e., volume of water densified from the upper to the lower limb of the AMOC) due to surface fluxes is associated with the SFWMT across the averaged AMOC isopycnal,   moc .
Maps of SFWMT are used to visualize the spatial distribution of the dense water formation due to surface fluxes over the subpolar gyre.These maps are estimated by integrating the buoyancy forcing given in Equation 1 over the outcropping area of specific isopycnals that bracket an isopycnal σ.For dense water formation due to surface fluxes, we consider the isopycnals that bracket   moc across each section with a density bin of 0.1 kg m −3 .The SFWMT is then divided by the grid area at each grid point (Sv m −2 ) to account for the difference in resolution between the models.
Finally, the dense water formation can be forced by interior mixing.This component, called Residual, is estimated as the difference between the overturning divergence over the Irminger and Iceland Basins and the SFWMT IIB across   moc (Evans et al., 2022).Over the Labrador Sea, the Residual is the difference between the overturning at OSNAP West and the SFWMT LS across   moc .

Evaluating the Simulated Overturning at OSNAP
The models are first evaluated by comparing the simulated overturning stream function with the OSNAP observations in terms of mean state and variability.
The time-averaged overturning at OSNAP East and OSNAP West in the models are broadly consistent with the observations (Figure 1).In particular, the majority of the overturning occurs at OSNAP East instead of OSNAP West in all models, which is consistent with the results of Menary et al. (2020) and Jackson and Petit (2022).However, significant biases remain.The magnitude of the overturning at OSNAP in the MM and HH models are 25% and 40% larger than in the observations.More precisely, the MM and HH models are both ∼4 Sv larger than the observed 13.9 ± 2.8 Sv at OSNAP East, and 2.3 and 5.3 Sv larger than the observed 3.1 ± 1.6 Sv at OSNAP West, respectively (Table 1).On the contrary, the LL model is 1 Sv smaller at OSNAP East and ∼2 Sv smaller at OSNAP West than in the observations.The net overturning in LL is thus statistically weaker than in the MM and HH models at the p < 0.05 level for the two sections (based on a Student's t-test), with a difference of 8.8 Sv between the HH and LL models at the full OSNAP line (Figure S2 in Supporting Information S1).

Note.
The overturning divergence between GSR and OSNAP East has been estimated as the maximum of the difference between the overturning stream functions of these two sections.We use  1), while the observed values are of 27.55 and 27.7 kg m −3 , respectively.The AMOC isopycnal is thus shallower along the OSNAP sections in the models than in the observations (not shown).
Because the OSNAP observations spans over 3 full years, we now evaluate the monthly variability of the overturning in the model simulations by comparing the moving standard deviation of the simulated overturning over periods of 3 years to the observed standard deviation.Figure 1c shows that the observed variability falls inside the distribution of the 3 models at OSNAP East, and inside the 25th to 75th percentiles for the HH model.However, the MM and HH models clearly overestimate the overturning variability at OSNAP West.The magnitudes of the overturning variability at OSNAP East and West are thus similar in these model simulations, which is not consistent with the OSNAP observations where OSNAP East dominates the variability of the full OSNAP section.Part of this difference is due to a large seasonal variability of the overturning at OSNAP West in the MM and HH models (0.3 Sv in LL, 1.7 Sv in MM, and 1.5 Sv in HH, not shown).Nevertheless, we note that the seasonal variability of the overturning at OSNAP in the 3 models is maximum in spring and minimum in winter, which is consistent with OSNAP observations (Li et al., 2021).

The Role of SFWMT in Explaining the Overturning Biases
Given the disagreements in overturning between models and observations, we now further investigate the causes of the biases in the models.In particular, we compare the overturning simulated at OSNAP to the dense water formation induced by surface fluxes north of the OSNAP arrays.From here on, we focus on the time mean overturning and SFWMT.
To locate the origin of the biases at OSNAP East, we first examine the strength of the overflow transport across the GSR.We find that the GSR stream function is similar between the 3 models with a maximum overturning ranging from 4.4 Sv in LL to 5.8 Sv in HH (Figure 1b), which is somewhat smaller than the observed overflow transport of 6.6 ± 0.4 Sv.Hence, the large overturning in MM and HH is not associated with differences in the overflow transport across the GSR and must be explained by larger dense water formation between the OSNAP line and the GSR.Such increased dense water formation can only be caused by increased surface fluxes or by interior mixing.
To estimate the overturning over the area bounded by OSNAP East and the GSR (i.e., the Irminger and Iceland Basins in Figure 2a), we compute the overturning divergence between the sections.The overturning divergence is more than 3 Sv larger in MM and HH than in LL, with notable differences in   moc between the 3 models (Table 1).In observations, the volume budget estimated between similar sections is of ∼7 ± 3.8 Sv in Petit et al. (2020), 9.6 ± 3.4 Sv in Chafik and Rossby (2019) and 10.2 ± 1.7 Sv in Sarafanov et al. (2012).Although the range of observed overturning is large due to the difference in observational data sets and time periods, we show that only the simulated overturning of the LL model falls inside this range.
Because the overturning is larger in MM and HH than in the observations or in LL, we now analyze the role of SFWMT for the overturning over the two areas (e.g., SFWMT IIB and SFWMT LS in Table 1).The SFWMT IIB peaks at a density of ∼27.5 kg m −3 over the Irminger and Iceland Basins (Figure 2b), which is associated with a densification of subpolar mode water over the Iceland Basin (Evans et al., 2022;Petit et al., 2020).Across   moc , the SFWMT IIB is of 8.2 ± 4.1 Sv in the LL model.This is close to the estimates from ECCO of 7.4 ± 1.6 Sv and from Petit et al. (2020) of 7.0 ± 2.5 Sv.However, the SFWMT IIB across   moc is larger in MM and HH with 10.4 ± 2.8 Sv-11.7 ± 2.2 Sv of dense water formed by surface fluxes in these models.The SFWMT IIB is thus mainly responsible for the large overturning found at medium and high resolutions over the Irminger and Iceland Basins.A similar picture emerges over the Labrador Sea in that the SFWMT LS explains most of the difference in overturning at OSNAP West between the models.Indeed, the SFWMT LS across   moc increases from 0.9 ± 1.0 Sv in LL to 5.1 ± 1.5 Sv in HH (Figure 2d), while estimates from ECCO and Petit et al. ( 2020) are of 2.0 ± 2.3 Sv and 1.5 ± 0.7 Sv, respectively.Moreover, water denser than 27.8 kg m −3 are further densified in the Labrador Sea for both the MM and HH models: the SFWMT LS peaks at 27.85 kg m −3 for MM (5.8 Sv) and HH (6.2 Sv).This densification of dense water is consistent with the simulated overturning stream function at OSNAP West that is positive at ∼27.85 kg m −3 for these models, while the observed stream function is negative at these densities (Figure 1e).
Finally, interior mixing possibly contributes to the increase in overturning with resolution, especially at OSNAP West.The contribution of the interior mixing to the total overturning is estimated as the residual between overturning and SFWMT across   moc (Table 1) and is found to be non negligible in the MM and HH models: up to ∼40% of the overturning over the Labrador Sea and ∼20% of the overturning over the Irminger and Iceland Basins.However, a comparison between the 3 models shows that the interior mixing clearly increases with resolution over the Labrador Sea, from 0.1 Sv in LL to 3.3 Sv in HH, while the residual in HH (1.4 Sv) is equal to that in LL and smaller than in MM (2.6 Sv) over the Irminger and Iceland Basins.Thus, in the model HadGEM3-GC3.1, the interior mixing should be less impacted by changes in resolution over the Irminger and Iceland Basins than over the Labrador Sea.

Linkage Between Local Biases and Intense SFWMT Over the Labrador Sea
Because the SFWMT is the leading contributor to the increasing overturning from low to high resolution, we now further analyze the local biases that influence the SFWMT.In particular, we focus on the SFWMT over the Labrador Sea because the overturning increases by ∼80% between the LL and MM models and by ∼40% between the MM and HH models across OSNAP West (p-value < 0.05).On the contrary, the overturning is similar between MM and HH over the Irminger and Iceland Basins, such that the mechanisms leading to this bias are not clearly affected by a difference in horizontal resolution from a ¼° to a 1/12° resolution.

Local Strengthening of SFWMT Over the Basin
The causes of the biases in SFWMT over the Labrador Sea are first investigated by decomposing the SFWMT LS into its thermal and haline components (i.e., each of the two terms in Equation 1).The thermal component of the SFWMT LS densifies surface water while the haline component decreases the density of surface water (Figure 3).Although the magnitude of the two fluxes increases with resolution over the Labrador Sea, the increase of the thermal surface density flux is relatively larger than that of the haline surface density flux.Therefore, the increase in SFWMT LS with resolution is dominated by a stronger increase of its thermal component that is not entirely compensated by its haline component.Furthermore, the haline and thermal components of the SFWMT LS are similarly distributed over the Labrador Sea, with an intensified strip of negative and positive transformations along the boundary of the basin (Figure 4).The fluxes are thus intensified along the sea ice edge in the Labrador Sea (as discussed in Section 5.2) such that the haline surface density flux is probably dominated by the melting/freezing of the sea ice.To quantify its impact for the SFWMT LS , the haline component of the surface density flux was computed without the melting/freezing term of the sea ice.Without this term, the haline surface density flux is close to zero at all densities for the 3 models (dashed blue line in Figure 3).Hence, the sea ice melting contributes significantly to the compensation of the increasing heat loss along the boundary of the Labrador Sea with resolution.
To better understand how the partial compensation of these two fluxes influences the total SFWMT over the Labrador Sea, we now compare the spatial pattern of SFWMT across   moc .In addition to a strengthening in SFWMT with increasing resolution, the location of dense water formation is shifted from the interior of the basin at low resolution to the boundary at higher resolutions (Figures 5b-5d).When using the 3,000-m isobath to bound the interior of the Labrador Sea, we estimate that 93.4% of SFWMT across   moc is localized over the interior of the basin in the LL model, while it is only 23.2% and 14.1% in the MM and HH models, respectively.
It is estimated to be 64.4% in ECCO (Figure 5a), which is consistent with the observed estimate in Li et al. (2021) from NCEP/NCAR and ERA5 surface fluxes.In contrast, a large densification of water at density classes denser than 27.8 kg m −3 occurs over the interior of the basin in the MM and HH models, whereas the SFWMT across these densities is negligible for ECCO (Figures 2d and 5e-5h).Therefore, the strengthening of SFWMT LS with resolution is associated with an increase of its thermal component along the boundary of the basin.

Role of Surface Density for Dense Water Formation
A shift in the location of dense water formation from the interior to the boundary of the basin can be related to a bias in density at the surface.Indeed, a modification in the spatial distribution of surface density influences the SFWMT by modifying the outcropping area of   moc over which the buoyancy flux is integrated (Equation 2).Petit et al. (2021) showed that this mechanism plays a key role for the variability in SFWMT over the eastern subpolar gyre, which is more sensitive to changes in the surface density field than to the direct influence of surface fluxes.
The surface distribution of   moc is compared between the models and ECCO alongside the thermal surface density flux in Figures 4a-4d.While the 27.65 kg m −3 isopycnal outcrops over the interior of the basin in ECCO and in the LL model, the isopycnal is greatly expanded over the boundaries of the Labrador basin in the MM and HH models.Similarly, while the denser 27.75 kg m −3 isopycnal does not outcrop over the Labrador Sea in ECCO, it does outcrop along the boundary of the basin in the MM and HH models.The outcropping area of   moc is thus shifted from the interior to the boundary of the basin with increasing resolution, such that the distribution of   moc in the MM and HH models matches that of intense surface density flux.Therefore, the shift in SFWMT across   moc over the Labrador Sea in these models is mainly attributed to anomalously dense surface water that shifts the outcropping area of   moc along the boundary of the basin (Figure 6 and Figure S1 in Supporting Information S1).The bias in surface density in the MM and HH models is not only found over the Labrador Sea, but also over the Irminger Sea (Figure 6).This bias is driven by a bias in salinity that brings too salty surface water from the subtropical gyre to the Irminger and Labrador seas (Jackson et al., 2020).This suggests that these biases do not only impact the dense water formation due to surface fluxes over the Labrador Sea (e.g., SFWMT LS across   moc ), but also its formation over the Irminger Sea.Furthermore, the profiles along OSNAP West show that these biases are found over the entire water column (Figure S3 in Supporting Information S1), which can modify the stratification and, thus, the internal mixing responsible for the dense water formation in the MM and HH models.
In contrast to the Irminger and Labrador seas, the surface salinity and density over the Iceland Basin in ECCO is close to that in the MM and HH models, while the LL model is too fresh over this basin.This suggests that the densification of subpolar mode water over the Iceland Basin is biased in the LL model instead of the MM and HH models.
Finally, a positive anomaly of surface temperature brings warm water over the Labrador Sea in the MM and HH models (Figure 6; Figure S3 in Supporting Information S1).In the following section, we discuss its impact on the representation of the sea ice and heat flux over the Labrador Sea.

Role of Sea Ice Extension for Dense Water Formation
The representation of the sea ice in the models impacts the dense water formation through different processes.First, the melting/freezing of the sea ice dominate the haline component of the surface density flux, as discussed in Section 5.1.However, sea ice also influences the dense water formation by insulating the ocean from the atmosphere and regulating air sea interactions during winter.In particular, the sea ice retreat during winter can lead to substantially enhanced turbulent heat loss and mechanical mixing (e.g., modification of the stratification via brine rejection) when warm boundary currents are exposed (Wu et al., 2021).Therefore, given the large surface density fluxes over the boundaries of the Labrador Sea (Figure 4), we address the question of how important the sea ice is for generating large SFWMT in the MM and HH models.

(e-h) To localize the transformation of water denser than
moc , the SFWMT has been integrated within the density bin 27.8-27.9kg m −3 for models and observations.SFWMT higher than 9 × 10 −12 is outlined in black.Gray line indicates the 3,000-m isobath from ETOPO1.Thick black line indicates the transect used for OSNAP West.The sea ice cover is analyzed alongside the heat loss over the Labrador Sea in Figure 7.In ECCO, the averaged Marginal Ice Zone (MIZ, i.e., area where the ocean has sea ice cover of 80%-15%) reaches the interior of the Labrador Sea in March (Figure 7a).The MIZ in the LL model appears broadly consistent with ECCO, although it extends further over the interior of the basin and covers the West Greenland Current region (Figure 7b).In contrast to LL and ECCO, there is substantially less sea ice in the MM and HH models.Indeed, the extension of the MIZ into the Labrador Sea is small in these models: the limit of 15% in sea ice concentration is localized much closer to the western boundary of the Labrador Sea during most of the 100 years of simulation, and the West Greenland current is largely ice free (Figures 7c and 7d; Figure S4 in Supporting Information S1).
Consequently, there appears to be a strong relationship between the distribution of sea ice concentration and heat loss over the Labrador Sea.To highlight this relationship, we compare their evolution along a section perpendicular to the edge of the sea ice (black line in Figures 7a-7d).In the MM and HH models, the heat loss clearly peaks at ∼60°W and is associated with an abrupt reduction in sea ice concentration (i.e., the edge of the sea ice) along the boundary of the basin, which is reinforced by a sea ice extent particularly constant over time in these models (Figures 7g and 7h; Figure S4 in Supporting Information S1).On the contrary, the larger MIZ in the LL model reduces the averaged heat loss that does not peak over the Labrador basin (Figure 7f), which is close to the heat loss distribution in ECCO (Figure 7e).Furthermore, the small MIZ area and intense heat loss in the MM and HH models are co-located with fast upper ocean velocities.Indeed, the comparison of meridional velocity with the sea ice concentration clearly shows that both models have a relatively intense boundary current at ∼60°W (Figures 7i-7l).The boundary current peaks at a similar longitude in the LL model, but is much weaker and broader in contrast.Thus, the large heat loss in MM and HH is consistent with an intense and warm boundary current (Figures 6e-6h) that is exposed to the atmosphere due to a lack of sea ice (which, in turn, is likely due to the warmer boundary currents).In contrast, the boundary currents in LL are cooler and insulated from the atmosphere by sea ice.

Discussion
In this study we have evaluated the subpolar AMOC in a hierarchy of HadGEM3-GC3.1 and have explored the reasons for biases in the simulated overturning at OSNAP.However, there are some remaining issues that warrant  further discussion, including the importance of sea ice for the overturning, the origin of the surface biases over the subpolar gyre, the role of interior mixing for the overturning, and the wider implications of these results as compared to other high-resolution models.
First, our study reveals that the representation of the sea ice plays a major role in the dense water formation over the Labrador Sea by impacting surface density fluxes along the boundary of the basin.Although our study only focuses on its impact in the mean state, Aylmer et al. (2022) showed that heat transport variability is a major contributor to the long-term evolution of sea surface temperature and sea ice extent over the Arctic in CMIP6 models.While they did not make a direct connection between heat transport variability and AMOC variability, their results and ours suggest a tight two-way coupling between ocean circulation and sea ice through SFWMT.Over the Labrador Sea, we also show that a bias in surface salinity is mainly responsible for the shift in dense water formation induced by the air-sea fluxes, intensified along the boundary of the basin in MM and HH.It raises the question of whether this shift led to a modification of the exchanges between the interior and boundary of the basin, which can possibly impacts the mechanisms of density compensation described by Zou et al. (2020) at OSNAP West.Moreover, a bias in salinity over the Labrador Sea led to a stronger AMOC weakening, as shown by Jackson et al. (2020).
Over the Irminger and Iceland Basins, it is still unclear if an increase of resolution from MM to HH influences the overturning over the area bounded by the GSR and OSNAP East sections.Because a similar density bias is found over the Irminger Sea (Figure 6), we hypothesize that the two models have a similar shift in SFWMT over the basin and that it possibly leads to a similar overturning divergence between the medium and high resolutions (Table 1).A deeper analysis of the impact of resolution on the representation of the overturning over the eastern subpolar gyre is thus needed.Overall, the origin of the biases in surface fields (e.g., temperature, salinity, and density) is difficult to assess because of the numerous interconnections between the biases in these variables as well as the complex interconnections between the subpolar and subtropical gyres (Kostov et al., 2021).The biases in surface salinity and temperature are possibly advected from the subtropical gyre by a misrepresentation of the large-scale circulation over the North Atlantic.Indeed, Jackson et al. (2020) suggested that anomalously salty water in high-resolution coupled models compared to low resolution models could be related to a stronger subpolar gyre and better position of the North Atlantic current.However, an improvement of the North Atlantic current position is not necessarily associated with an improvement of the subpolar gyre circulation.Indeed, Treguier et al. (2005) showed that models with a ¼° ocean resolution had a good position of the North Atlantic current but too much transport across the Reykjanes Ridge and too little freshwater export in the boundary currents from the Arctic.The smaller bias in the overturning found across OSNAP in LL compared to the higher resolution models is, thus, not necessarily associated with a better representation of the North Atlantic, but is rather related to a compensation of biases advected from the subtropical gyre.
One caveat to consider in the study is the short spin-up of the 3 models, and the associated drifts, that can be a source of uncertainty for the intercomparison of the models.Rattan et al. (2010) showed that the drift in salinity over the Labrador Sea rapidly developed during the first 3-years period of the initial adjustment and is mainly due to the representation of freshwater in the boundary currents and to the advection of salty subpolar mode water.In line with this study, the salinity and temperature conditions during the 30 years spin-up simulation of our 3 models shows that salty and warm surface water quickly spreads over the subpolar gyre in the MM and HH models from the subtropical gyre (Figure S5 in Supporting Information S1).The decrease in density over the subpolar gyre during the spin-up of the LL model is also consistent with previous studies showing a northward advection of fresh and cold anomalies in low resolution models (Athanasiadis et al., 2022;Roberts et al., 2020).Nonetheless, the magnitude of the drifts in thermohaline properties caused by a short spin-up is considerably smaller than the biases observed between the models, which is consistent with Roberts et al. (2019).Furthermore, the simulated low frequency variability in overturning at OSNAP, estimated as the standard deviation of the moving mean in overturning during periods of 50 years (0.9 Sv for LL; 0.2 Sv for MM; 0.4 Sv for HH) is also small compared to the 8.8 Sv of difference between the 100-years mean states of the high and low resolution (Figure S2 in Supporting Information S1).While long term drifts may develop in a non-linear way, it appears that the short spin-up in our simulations only explains a small fraction of the differences between the modeled AMOC mean states.
Moreover, although we focused on SFWMT, other processes related to dense water formation over the subpolar gyre can be biased in these models; such as interior mixing.One way to assess the interior mixing in the models is by examining the thermal and haline tendencies due to diffusive fluxes, which were saved for the LL and MM models (Figures S6 and S7 in Supporting Information S1) but not for the HH model.The overall distribution of these fluxes is similar between the two models, intensified along the boundary of the Irminger and Labrador seas at sub-surface.However, the fluxes over the Labrador Sea are stronger along the boundary of the basin and are ∼10 times stronger above 500-m depth in MM than in LL (Figures 7j and 7k).The difference in diffusive fluxes can be related to differences in vertical velocity.Indeed, Katsman et al. (2018) hypothesized that a larger diffusivity of a 1° model can allow sinking to occur in the interior of the basin by breaking geostrophy as compared to a higher resolution model.
Finally, our study shows that the overturning at OSNAP increases with resolution in the HadGEM3-GC3.1 model and is too large at medium and high resolutions as compared to observations.A key question is whether these biases have a wider effect on the North Atlantic.In particular, to what extent the subpolar overturning biases in MM and HH contribute to the larger shift in the Gulf Stream separation observed by Grist et al. (2021) in eddy-rich simulations, or to the larger AMOC forced response at 26.5°N in the higher-resolution models found by Roberts et al. (2020).Using the same hierarchy of HadGEM3-GC3.1,Roberts et al. (2019) found a similar increasing overturning with increasing resolution at 26.5°N.The simulated overturning in HH was stronger than observations at the RAPID array.However, the overturning was close between LL and MM and weaker than observations.In contrast, our analysis of the subpolar AMOC revealed similar overturning biases between the MM and HH models at the OSNAP line.This shows a non-trivial link between AMOC biases at different latitudes, suggestive of different underlying causes.Furthermore, the resolution dependence of the subpolar overturning and SFWMT in HadGEM3-GC3.1 differs from that seen in the CESM1 model, where the simulated overturning and SFWMT is much reduced and closer to observations at high resolution as compared to the low resolution, especially at OSNAP West (Oldenburg PETIT ET AL.  et al., 2022;Yeager et al., 2021).Hence, although models at high resolution better represent small-scale processes (e.g., eddies, fronts, boundary currents…) than model at low resolution (Fox-Kemper et al., 2019), the analysis presented here shows that increased resolution does not lead to a systematic improvement of the subpolar overturning by itself.Rather, our study has underlined a broad range of interactions that can affect its simulation at high resolution.The associated causes can be diverse (e.g., spurious mixing on the grid scale or the need for finer resolution to fully resolve the Rossby radius at these latitudes) and needs to be further analyzed.It should also be noted that, although the HighResMIP protocol has been designed to evaluate the sensitivity of the simulated global climate to horizontal resolution, we cannot attribute changes in ocean processes only to changes in ocean resolution because of the strong interconnection between ocean, ice and atmosphere.Thus, further experiments would be required to evaluate the direct impact of resolution on the biases identified here.

Conclusions
In this study, a resolution hierarchy of the global coupled model HadGEM3-GC3.1 was used to evaluate the sensitivity of dense water formation to horizontal resolution over the subpolar gyre.We first evaluated the overturning of 3 models with 100/250 km (LL), 25/100 km (MM), and 8/50 km (HH) ocean/atmosphere resolutions in comparison to OSNAP observations.The part of dense water formation due to surface fluxes (i.e., the surface forced water mass transformation, SFWMT) was then estimated north of the OSNAP line in the models and compared to the ocean reanalysis ECCO.The main results are summarized below.
The simulated overturning stream functions at OSNAP are broadly consistent with the observations in all resolutions as most of the overturning occurs at OSNAP East instead of OSNAP West.However, the strength of the overturning at OSNAP is significantly larger in the MM and HH resolutions (25% and 40%, respectively) than in the observations.The increase in overturning with resolution is particularly clear at OSNAP West, where the overturning in the MM and HH is 170% and 270% larger than in the observations, respectively.Furthermore, the monthly variability of the overturning at OSNAP West is equal or more variable than at OSNAP East in the MM and HH models, which is inconsistent with the OSNAP observations where OSNAP East dominates the variability of the full OSNAP section.
The disagreement between these models and observations is not attributed to anomalously large overflow transport through the GSR.Rather, the disagreement is consistent with anomalously large SFWMT over both the Labrador Sea and the Irminger and Iceland Basins.Interior mixing in the models at high resolution also contributes to up to ∼40% and ∼20% of the total overturning over these two areas, respectively.Over the Labrador Sea, the dense water formation due to both surface fluxes and interior mixing increases between the three resolutions.Over the Irminger and Iceland Basins, the dense water formation increases only between the LL and MM resolutions but not between the MM and HH resolutions.
North of OSNAP West, the intensification in SFWMT LS at   moc with resolution is attributed to a combination of two main biases over the Labrador Sea: a bias in surface temperature and a bias in surface salinity that are advected by the boundary current from the eastern subpolar gyre.The warming is responsible for reducing the marginal sea ice along the boundary of the basin and reducing the area of insulation of the ocean, which enhances heat loss from the exposed boundary current.Additionally, anomalously salty surface water leads to anomalously dense surface water over the Labrador Sea.This bias shifts the outcropping area of   moc from the interior to the boundary of the basin, where the largest surface density fluxes occur, and thus results in intense SFWMT over the Labrador Sea at MM and HH resolutions.Therefore, our analysis shed light on a range of model biases related to an increase in horizontal resolution of the model HadGEM3-GC3.1 that needs to be improved in order to better simulate the overturning over the Labrador Sea.However, it is unclear how these biases have a wider impact on the overturning further south or on the AMOC variability and, in particular, on the forced response of the AMOC investigated at 26.5°N by previous works.The origin of these biases in the subpolar gyre and their impacts on the overturning over the eastern subpolar gyre are also unclear and needs to be further analyzed.

Figure 1 .
Figure 1.Comparison of overturning stream function between models and observations.(a) Map of localization of the sections: OSNAP West, OSNAP East and Greenland-Scotland Ridge (GSR).The overturning stream functions are shown for (b) GSR, (d) the entire OSNAP line, (e) OSNAP West and (d) OSNAP East.The black lines are observations from the OSNAP array.The black circle and bar in panel (b) indicates the averaged overflow transport and its uncertainty at the GSR of 6.6 ± 0.4 Sv.Note that the vertical positioning of the transport at a density of 27.8 kg m −3 is arbitrary and associated with the upper limit of the overflow layer.Shaded areas indicate their interannual standard deviations.(c) Box plot of the moving standard deviation of the overturning at OSNAP West and East.The moving standard deviation is estimated from the monthly values at the maximum overturning stream functions over 3-years windows and includes the seasonal variability.Similarly, the black crosses indicate the standard deviations of the monthly overturning over the 3 full years of OSNAP observations, of 2.9 Sv at OSNAP East and 1.4 Sv at OSNAP West.The red line indicates the median and the edges of the blue box indicate the 25th and 75th percentiles.The whiskers indicate values within 1.5 times the interquartile range, while the red crosses indicate outlier values that are outside of this range.

Figure 2 .
Figure 2. Overturning and SFWMT stream functions over the subpolar gyre.Over the Irminger and Iceland Basins, we compare the (a) overturning divergence between OSNAP East and GSR with the (b) SFWMT IIB .Over the Labrador Sea, we compare the (c) overturning stream function at OSNAP West (as shown in Figure 1e) with the (d) SFWMT LS .The SFWMT estimated with ECCO are indicated by black lines.Shaded areas indicate the interannual standard deviation.

Figure 3 .
Figure 3. Decomposition of the SFWMT LS into its (red) thermal and (blue) haline contributions for the (a) LL, (b) MM, and (c) HH models.The dashed blue lines show the haline contribution of the SFWMT LS estimated without the melting/freezing term of the sea ice.

Figure 4 .
Figure 4. (a-d) Thermal and (e-h) haline contribution of the surface density flux (10 −11 Sv m −2 ).The thin dashed and plain lines in the top panels indicate the outcropped isopycnals 27.65 and 27.75 kg m −3 in March, respectively.Gray line indicates the 3,000-m isobath from ETOPO1.Thick black line indicates the transect used for OSNAP West.

Figure 5 .
Figure 5. Spatial distribution of the SFWMT estimated across specific isopycnals and divided by the area of each grid (10 −12 Sv m −2 ).To localize the dense water formation due to surface fluxes, the SFWMT has been integrated within the density bin associated with   moc at OSNAP West: (a) 27.65-27.75kg m −3 for observations and (b-d) 27.7-27.8kg m −3 for models.(e-h)To localize the transformation of water denser than  moc , the SFWMT has been integrated within the density bin 27.8-27.9kg m −3 for models and observations.SFWMT higher than 9 × 10 −12 is outlined in black.Gray line indicates the 3,000-m isobath from ETOPO1.Thick black line indicates the transect used for OSNAP West.

Figure 7 .
Figure 7. Linkage between heat loss and marginal ice zone.(a-d) Heat loss averaged in March (W m −2 ).The black plain (80% of ice concentration) and dashed (15% of ice concentration) lines indicate the Marginal Ice Zone averaged in March.(e-h) Surface heat loss (blue lines) and (i-l) Meridional velocity averaged in the upper 100-m depth (blue lines) as compared with Sea Ice Concentration (SIC; black lines) integrated in March along the section at 58°N indicated by the black lines in panels (a-d).The black shaded areas indicate the standard deviation of the SIC.