Figure 1A presents the time-mean (2009–2018) ocean surface stress in CONTROL, which is characterized by large values in the western boundary region of the Labrador Sea (45°N–65°N, 30°W–45°W), the Greenland Sea (65°N–80°N, 45°W–0°), and the North Atlantic Ocean (50°N–65°N, 45°W–0°); small values in the central region of the Nordic Seas (65°N–80°N, 45°W–15°E). Its pattern is similar to results from the previous numerical simulations (Martin et al., 2014; Wu et al., 2016). When the ocean surface currents are excluded from the ice-ocean stress calculation, an increase of time-mean ocean surface stress occurs in regions covered by sea ice, especially in the Baffin Bay, the western boundary areas of the Nordic Seas, and the Labrador Sea (Figure 1B). In the ice-free region, the non-zero values are induced by the different ocean surface currents, which are used to calculate the air-sea stress (Figure 1B). This increase in ocean surface stress is more pronounced during the winter (not shown). Quantitatively, the magnitude of the time-mean ocean surface stress averaged over the whole (sea ice-covered) region is about 7.7% (18%) weaker in CONTROL than that in NONE (Table 1). Additionally, a reduction of 11 and 15% in ocean surface stress was also found in the Labrador Sea and the Nordic Seas, respectively. The significant increase in the time-mean ocean surface stress indicates that sea ice drift is indeed generally orientated in directions similar to the ocean surface currents (Figure 1B).

As shown in Figure 1D, the spatial pattern of Ekman pumping is characterized by large values in the Nordic Seas, the Irminger Sea, and the Labrador Sea, especially in the region around the western boundary of the Greenland Sea (Figure 1D). Ekman pumping was calculated using k⋅∇×(τρf), where k is the unit vector in the vertical direction, τ is the ocean surface stress, f is the Coriolis parameter, and ρ is the seawater density. The increase in the time-mean ocean surface stress also leads to a widespread enhancement in the magnitude of Ekman pumping in NONE, most noticeable in the Baffin Bay, the Labrador Sea, and the western boundary region of the Greenland Sea (Figure 1E). For example, the strength of the time-mean Ekman pumping averaged over the Nordic Seas is about 19% stronger in NONE than that in CONTROL; moreover, an increase of about 17% is also found in the Labrador Sea (Figure 1E).

Due to the scarcity of observations, the ocean surface currents were often replaced by surface geostrophic currents in previous studies when calculating ice-ocean stress (Yang, 2006, 2009; Zhong et al., 2018). Here, the different effects of including total ocean surface currents (ageostrophic and geostrophic components) or only including geostrophic currents on ice-ocean stress and Ekman pumping have also been examined (Figures 1C,F). Reducing the surface currents to their geostrophic component significantly impacts ice-ocean stress, especially in the Nordic Seas and the Labrador Sea (Figures 1B,C). The pattern and magnitude of ice-ocean stress calculated by the geostrophic component differ from that calculated by the total ocean surface currents, especially in the Labrador Sea, the Irminger Sea, and the Baffin Bay where the difference virtually disappears (Figures 1B,C). Similar results are also found in the Ekman pumping (Figures 1E,F). These above results indicate that the ageostrophic currents play an important role in the calculation of ice-ocean stress and Ekman pumping. When calculating the ice-ocean stress, the ageostrophic current should be properly considered in further studies rather than only the geostrophic current as done in previous studies (Yang, 2006, 2009; Zhong et al., 2018).

Figure 2 shows the time-mean net surface freshwater flux (the sum of evaporation, precipitation, runoff, and sea ice freezing/melting), net heat flux (the sum of sensible heat flux, latent heat flux, long-wave radiation, and short-wave radiation), SST, and sea surface salinity (SSS) averaged over the last decade in CONTROL and their differences between CONTROL and NONE. The overall patterns of time-mean surface freshwater and heat fluxes are very similar in both experiments, with an extensive heat loss in the North Atlantic and a significant freshwater input in the western boundary regions of the Greenland Sea and the Labrador Sea, seasonally covered by the sea ice (Figures 2A,D). In sea ice-covered regions, there is a significant enhancement in surface heat loss and freshwater input when ocean surface currents are excluded from ice-ocean stress calculation (Figures 2B,E). Heat loss (freshwater input) averaged over the sea ice-covered region in Figure 2B decreases by 54.4% (75%) from −79 W m^–2 (8 × 10^–4 kg m^–2 s^–1) in NONE to −36 W m^–2 (2 × 10^–4 kg m^–2 s^–1) in CONTROL (Table 1). However, the heat loss (freshwater input) in the North Atlantic Ocean in CONTROL decreased by 12.1% (54.5%) compared with that in NONE (Table 1). The differences in surface freshwater and heat fluxes between CONTROL and NONE are closely linked to their differences in SST and SSS (Figures 2C,F). The significant heat flux differences in the North Atlantic are in opposition to the SST differences, which demonstrates that changes in SST lead to changes in surface heat flux as well (Figures 2B,C). The following will show that the SST differences between CONTROL and NONE in the North Atlantic Ocean are mainly associated with differences in the strength of northward heat transport. As shown in Figure 2F, the SSS increases significantly over the sea ice-covered regions as a result of the stronger upwelling in NONE than that in CONTROL (Figures 1B,E), even with more freshwater entering into the ocean (Figure 2E). The freshwater flux induced by the sea ice melting\freezing is almost the same as the net freshwater flux (not shown). This confirms that more freshwater input in NONE comes from increased SST-induced sea ice melting than that in CONTROL. However, a slight decrease in SSS is also found in North Atlantic Ocean, which is caused by the spread of freshwater in NONE (Figures 2E,F and Table 1).

Figure 3A shows the time-mean distributions of the mechanical energy input to the time-mean currents by the time-mean ocean surface stress (P). It is calculated using P=τ¯⋅u¯, where the over bar denotes a 10-year mean, τ is the ocean surface stress, and u is the ocean surface currents. The spatial patterns of P in both experiments are similar to those in previous studies (von Storch et al., 2012; Wu et al., 2017b), with most of the mechanical energy input concentrated in the central North Atlantic Ocean and the western boundary regions of the Greenland Sea, Labrador Sea, and Irminger Sea (Figure 3A). In NONE, the increase in P is significant, especially in the western boundary regions of the Irminger Sea, Labrador Sea, and Greenland Sea (Figure 3B). This increase is much more pronounced in wintertime (not shown). Figure 3C shows the distribution of mechanical energy input to time-varying ocean velocity by time-varying ocean surface stress (P′=τ′⋅u′¯) in CONTROL, which is an energy source of eddy kinetic energy (EKE) and characterized by large values in the Denmark Strait, the North Atlantic Ocean, and the Labrador Sea. Similarly, the difference in P′ between CONTROL and NONE is most significant in the western boundary regions of the Irminger Sea, the Labrador Sea, and the Greenland Sea, where there are the seasonally sea ice-covered regions (Figure 3D).

The spatial integrated value of P(P′) decreases by 39.0% (41.7%) from 0.41 TW (1 TW = 10¹² W; 0.24 TW) in NONE to 0.25 TW (0.14 TW) in CONTROL. Considering the ocean surface currents in the ice-ocean stress calculation results in a decrease in P (P′) by 0.16 TW (0.1 TW), of which 0.1 TW (0.07 TW) is due to a decrease in τ¯⋅u¯ag (τ′⋅uag′¯), and only 0.06 TW (0.03 TW) is due to a decrease in τ¯⋅u¯g (τ′⋅ug′¯), where the prime denotes the derivation from the 10-year mean; u_ag and u_g are the ageostrophic and geostrophic components of the ocean surface currents, respectively.

The effects of considering ocean surface currents in the ice-ocean stress calculation on the EKE and mean kinetic energy (MKE) are also investigated. The EKE and MKE are defined as EKE=0.5⋅(u′+v′)¯ and MKE=0.5⋅(u¯2+v¯2)’ where u and v are the zonal and meridional ocean currents, respectively. Similar to the previous studies (Eden and Böning, 2002; Rieck et al., 2019), the distributions of the time-mean MKE and EKE are characterized by large values along the boundary currents over the Irminger Sea, the Labrador Sea, the Greenland Sea, and in the vicinity of the North Atlantic Current branches (Figures 4A,C). When ocean surface currents are excluded from the ice-ocean stress calculation, there is a resultant increase in the surface MKE over the North Atlantic Ocean, most notably in the western boundary regions of the Labrador Sea, the Irminger Sea, and the Greenland Sea (Figure 4B). There is an EKE increase in the western boundary region of the Greenland Sea and the Baffin Bay where there are seasonally sea ice-covered regions (Figure 4D). It will be shown later that this EKE increase is due to the combined effect of the energy conversion from eddy available potential energy (EAPE) and MKE to EKE, which is induced by a stronger mechanical energy input (Figures 1B,E, 3B,D). However, the significant reduction of EKE is found in NONE across the North Atlantic Ocean and Nordic Seas (Figure 4D), which is due to the stable water column induced by more freshwater input in NONE than that in CONTROL, as previously shown (Zhong et al., 2018). Quantitatively, integrated MKE over the entire region decreases by about 33.3% from 0.3 EJ (1 EJ = 10¹⁸ J) in NONE to 0.2 EJ in CONTROL (Table 1). Spatially integrated over sea ice-covered regions, EKE in CONTROL and NONE are 4.5 PJ (1 PJ = 10¹⁵ J) and 8 PJ, respectively, representing a 43.8% reduction from CONTROL (NONE minus CONTROL) (Figure 4D). Although the magnitude of the MKE and EKE reduction weakens with depth, the reduction percentage of the area-averaged MKE and EKE remains roughly 10 and 25% deeper than 300 m (not shown).

The energy budget of EKE and MKE has also been calculated to explain the change in kinetic energy when considering the ocean surface currents in the ice-ocean stress calculation. The formulas used to calculate the energy conversion between EKE and MKE are the same as those used in Wu et al. (2017a). Figure 5A presents the distribution of vertically integrated energy conversion from EAPE to EKE through the baroclinic pathway. The energy is converted from EAPE to EKE (positive value) in the North Atlantic Ocean, the Labrador Sea, and the Nordic Seas (Figure 5A), which is similar to the previous studies (Eden and Böning, 2002; von Storch et al., 2012; Rieck et al., 2019). In NONE, a strengthened Ekman pumping leads to a more titled isopycnal that contains more available potential energy (APE) in the western boundary regions of the Nordic Seas and the Baffin Bay, which is released through the baroclinic pathway (Figure 5B). However, there is a reduction in the energy conversion from EAPE to EKE in the broad North Atlantic Ocean and the Labrador Sea (Figure 5B). Enhanced sea ice melting results in more freshwater input to the Labrador Sea and the North Atlantic Ocean in NONE than in CONTROL (Figures 2E,F), so stratification is much more stable there. Spatially integrated, the energy conversion from EAPE to EKE is 20 GW (1 GW = 10⁹ TW) in CONTROL, while it decreases to 12 GW in NONE, representing a 66.7% reduction.

In addition, the EKE can also obtain energy directly from the MKE through the barotropic pathway. Figure 5C gives the vertically integrated distribution of energy conversion from MKE to EKE, with a large energy conversion rate from MKE to EKE (positive values) in the Labrador Sea, especially in the west Greenland current region, which is similar to the previous studies (Eden and Böning, 2002; von Storch et al., 2012; Wu et al., 2016, 2017b; Rieck et al., 2019). When ocean surface currents are excluded from the ice-ocean stress calculation, the shear instability between sea ice drift and ocean surface velocity is strengthened in the sea ice-covered regions, such as in the western boundary regions of the Nordic Seas, the Irminger Sea, and the Labrador Sea (Figure 5D). Integrated over the North Atlantic Ocean (sea ice-covered region), the energy conversion from MKE to EKE is 1.9 GW (2.1 GW) in CONTROL, while it decreases to 1.1 GW (2.3 GW) in NONE, accounting for a 72.7% (−8.7%) decrease.

The changed air-sea fluxes and SST induced by the increased ice-ocean stress in NONE can be fed back to the overlying sea ice. As shown in Figure 6A, the SIC decreases extensively in NONE, especially in the western boundary regions of the Baffin Bay, the Labrador Sea, the Irminger Sea, and the Greenland Sea (Figure 6A). Quantitatively, the annual mean SICs in CONTROL and NONE are 33 and 22%, respectively, representing a 33.3% decrease with respect to CONTROL; the winter mean values are 47 and 22%, respectively, representing a 53.2% decrease with respect to CONTROL. In NONE, the sea ice extent (SIE) also expands dramatically, as shown by the labeled lines of 15% SIC in Figure 6, especially in the Nordic Seas (Figures 6A,B). Similarly, the sea ice thickness (SIT) also decreases significantly as large as 0.2 m in the western boundary regions of the Nordic Seas (Figure 6C); this decrease is much larger in the wintertime, amounting to 0.5 m (Figure 6D). Besides, there is a time-mean reduction of 27.6% (0.76 m in CONTROL and 0.55 m in NONE) and 29.6% (0.81 m in CONTROL and 0.57 m in NONE) of SIT from CONTROL to NONE throughout the whole year and wintertime, respectively (Figures 6C,D).

Figures 7A,B show the monthly mean total sea ice volume (SIV) and SIE in both simulations over the Arctic Ocean (regions where the SIC is equal to or greater than 15%). The SIE is much larger in NONE than in CONTROL. Quantitatively, the time-mean SIE in September decreases from 5.5 × 10⁶ km² in NONE to 4.5 × 10⁶ km² in CONTROL over the Arctic Ocean. As analyzed above, the significant increase in SIE is attributed to increased sea ice transport caused by greater ice-ocean stress in NONE. However, the total monthly SIV over the Arctic Ocean decreases entirely in NONE (Figure 7B). For example, the SIV in September decreases from 1.9 × 10⁴ km³ in CONTROL to 1.8 × 10⁴ km³ in NONE. The reduction in SIV is caused by a much thinner SIT in NONE than in CONTROL (Figures 6C,D). As shown above, the changes in SIE and SIV are the combined effect of increased ice-ocean stress and higher SST.

The mean kinetic energy (SIMKE) and eddy kinetic energy (SIEKE) of sea ice over the Arctic Ocean are also presented in both simulations (Figures 7C,D). The SIMKE and SIEKE are defined as SIMKE=0.5⋅(u¯i2+v¯i2) and SIEKEi=0.5⋅(ui′+vi′)¯, where u_i and v_i are the zonal and meridional sea ice velocity, respectively. In NONE, the movement of the sea ice becomes more active than that in CONTROL, and both SIMKE and SIEKE increase visibly. For example, the time-mean and area-mean SIMKE decreases from 0.0047 m^–2 s^–2 in NONE to 0.0039 m^–2 s^–2 in CONTROL (Figure 7C); besides, SIEKE also decreases greatly from 0.01 m^–2 s^–2 in NONE to 0.007 m^–2 s^–2 in CONTROL (Figure 7D).

Differences in the depth-integrated meridional volume transport and the gyre circulations between NONE and CONTROL are likely to be induced by the differences in the time-mean Ekman pumping. Figure 8A presents the time-mean barotropic stream function in CONTROL, which is similar to previous studies, characterized by cyclonic gyres in the North Atlantic Ocean and the Nordic Seas (Wu et al., 2016, 2017b). Figure 8B shows that the strength of the simulated subpolar gyre reduces almost everywhere when considering ocean surface currents in the ice-ocean stress calculation. For example, the strength of the time-mean subpolar gyre in the North Atlantic Ocean decreases by about 18.1% from 45.3 Sv in NONE to 37.1 Sv in CONTROL, and that of the gyre in the Nordic Seas reduces by about 12.6% from 15.1 Sv in NONE to 13.2 Sv in CONTROL.

We now investigate the impacts of considering ocean surface currents in ice-ocean stress calculation on MLD in the Nordic Seas and North Atlantic Ocean. The MLD is defined similarly to previous studies (Downes et al., 2009; Liu et al., 2017; Wu et al., 2020) as the depth at which the potential density is 0.03 kg m^–3 larger than that at the surface. The time-mean patterns of the MLD and the differences between NONE and CONTROL are shown in Figures 9A–C. The simulated MLDs in both experiments are represented by deep MLD in the Labrador Sea and Nordic Seas (Figure 9A), which is similar to the results from previous investigations (Kara and Rochford, 2003; Holdsworth and Myers, 2015; Wu et al., 2016). It is noted that the MLD is deeper than the observations, which is induced by a higher surface salinity than that from the observations. Additionally, the model resolution is relatively coarse, especially at the subpolar latitudes, and thus likely underestimates the role played by eddies in the restratification process after deep convection events. Excluding ocean surface currents from the ice-ocean stress calculation, the MLD increases considerably in the Labrador Sea and the Nordic Seas (Figure 9B). This increase is also more significant in wintertime (Figure 9C). For example, the time-mean MLDs averaged in the Labrador Sea (Nordic Seas) are 2,750 and 3,450 m (210 and 407 m) in CONTROL and NONE, respectively, representing a 20.3% (48.4%) increase. This difference amounts to 34 and 61% averaged over the Labrador Sea and the Nordic Seas in wintertime (Figure 9C). As shown in Figure 1E, in CONTROL, considering ocean surface currents in the ice-ocean stress calculation reduces the ocean surface stress and wind stress curl in the subpolar North Atlantic. For example, Ekman pumping is about 12.1% weaker in CONTROL than that in NONE with respect to the average for the Labrador Sea. The reduced wind stress curl weakens the strength of Ekman pumping and results in less doming of the density surface and a weaker cyclonic circulation over the Labrador Sea in CONTROL than that in NONE (Figure 8B). This reduced effect of pre-conditioning then makes it harder for sea surface heat loss to overcome stratification and trigger deep-reaching convection in CONTROL (Marshall and Schott, 1999).

The changes in the intensity of deep convection between these two experiments in the Labrador Sea and the Nordic Seas correlate with the differences in the strength of AMOC and MOC in the Nordic Seas, as found in previous studies (Eden and Willebrand, 2001; Zhai et al., 2011; Zhai and Johnson, 2014; Wu et al., 2016). The definition of MOC is similar to that in Wu et al. (2016, 2020) by integrating the meridional velocity across the ocean basin from the western boundary (x_W) to eastern boundary (x_E), and from the ocean bottom at z = −h upward: ψ(y,z,t)=∫-hz∫xW(y,z)xE(y,z)v(x,y,z,t)dxdz. The maximum value of ψ is defined as the strength of the MOC. The structure of AMOC is familiar in both simulations, with light water transports northward and dense water transports southward (Figure 10A; Wunsch and Heimbach, 2013; Wu et al., 2016). As shown in Figure 10B, when ocean surface currents are excluded, the strength of the simulated AMOC is coherently enhanced at all latitudes by as much as 3 Sv, representing a 14.3% reduction (Table 1). Besides, Figure 11A shows the distribution of the MOC in the Nordic Seas, which shows a similar structure in both simulations, with the light water transports northward and the dense water transports southward. In NONE, the strength of the simulated MOC in Nordic Seas is coherently enhanced at all latitudes (Figure 11B). The maximum MOC strength in Nordic Seas decreases by about 31.0% from 4.2 Sv in NONE to 2.9 Sv in CONTROL. One of the roles of the ocean in the climate system is the transport of heat to high latitudes. The meridional heat transport (MHT) can be adequately approximated (Warren, 1999; Zhang et al., 2002; Msadek et al., 2013) by ρCp∫0L∫-H0vTdzdx, where C_p is the specific heat capacity of seawater, T is the potential temperature, H and L are the depth and width of the ocean, respectively. In the Atlantic Ocean, the structures of MHT in both simulations are similar to the observed and simulated results (not shown, Ganachaud and Wunsch, 2003; Wu et al., 2016). For example, the maximum MHT is found to be 0.93 PW (1 PW = 10¹⁵ W) in CONTROL at 25°N, which is comparable to the 1.07 ± 0.26 PW and 1.27 ± 0.15 PW estimated in the previous studies (Macdonald, 1998; Ganachaud and Wunsch, 2003). In comparison, the maximum MHT increases to 1.05 PW in NONE, representing a 13% enhancement. It is noted that the MHT in NONE is closer to observations than that in CONTROL. It is likely caused by the weaker AMOC in ECCO2 than that in observations (Ganachaud and Wunsch, 2003), as the strength of AMOC dominates the MHT in the Atlantic Ocean.