North Atlantic subsurface temperature response controlled by effective freshwater input in “Heinrich” events

a Key Laboratory of Meteorological Disaster of Ministry of Education, Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disaster, Nanjing University of Information Science and Technology, Nanjing, China b Atmospheric Science Program, Department of Geography, Ohio State University, Columbus, OH, USA c Open Studio for Ocean-Climate-Isotope Modeling, Pilot National Laboratory for Marine Science and Technology (Qingdao), Qingdao, China d Department of Earth and Environmental Sciences, University of Michigan, Ann Arbor, MI, USA e Computational Physics and Methods (CCS-2) and Center for Nonlinear Studies (CNLS), Los Alamos National Laboratory, Los Alamos, NM, USA f Physical Oceanography Laboratory, Ocean University of China, Qingdao, China g National Center for Atmospheric Research, Climate and Global Dynamics Laboratory, Boulder, CO, USA h Department of Atmospheric and Oceanic Sciences, Peking University, Beijing, China i College of Computer, National University of Defense Technology, Changsha, China

The northern North Atlantic (NNA) subsurface temperature in response to the slowdown of the Atlantic meridional overturning circulation (AMOC) is crucial for ice sheet calving and recovery of the AMOC in Heinrich events. Paleoclimate proxies and modeling studies suggest that the NNA subsurface exhibits a robust warming during Heinrich 1, but with a less clear response during Younger Dryas (YD). The mechanism for the potentially different subsurface responses has remained not well understood. Previous studies show different NNA subsurface response depending on hosing at different locations. Here, by examining a suite of "water-hosing" experiments with different hosing regions and intensities, we show that, regardless of the hosing location, NNA subsurface temperature response is determined by the effective freshening over the NNA deep convective regions through the competition between the warming associated with the suppressed vertical mixing and the cooling associated with the weakened AMOC heat transport. A weak effective freshening favors advective effect and, in turn, cooling, while a strong effective freshening strengthens the mixing effect and leads to a warming. Our results suggest that a cooling may occur in the NNA subsurface during YD.

Introduction
Ocean sediments have recorded abrupt climate changes known as Heinrich events during the glacial-deglaciation cycles, notably the Heinrich 1 (H1) and Younger Dryas (YD) during the last deglaciation (Heinrich, 1988). These Heinrich events have been suggested to involve the Atlantic meridional overturning circulation (AMOC) and meltwater discharges over the North Atlantic (Marcott et al., 2011;Barker et al., 2015;McManus et al., 2004). A meltwater discharge reduces AMOC and the associated northward ocean heat transport, generating a robust surface bipolar seesaw response over the globe that is characterized by cooling in the Northern Hemisphere and warming in the Southern Hemisphere (Stocker and Johnsen, 2003;Stenni et al., 2011;Knutti et al., 2004;Buizert et al., 2015). In the subsurface (300 m to 2000 m) however, ocean temperature response in H1 tends to show a robust warming from the Southern Ocean into the tropical Atlantic and further into the northern North Atlantic (NNA), in both reconstructions (Marcott et al., 2011;Pedro et al., 2018;Weldeab et al., 2016;Barker and Diz, 2014;Dokken and Jansen, 1999;Rasmussen and Thomsen, 2004) and climate models (Liu et al., 2009;Zhu et al., 2017;Mignot et al., 2007;Brady and Otto-Bliesner, 2011;. While most proxy records and model simulations seem to infer a NNA subsurface warming in H1 compared to the last glacial maximum (Marcott et al., 2011;Thiagarajan et al., 2014;Petersen et al., 2013;Zhang et al., 2017;Liu et al., 2009), paleoclimate records seem to show some inconsistences in NNA subsurface tempera- ture response at YD, with some suggesting warming (Ezat et al., 2014) and some others cooling (Thiagarajan et al., 2014;Skinner et al., 2003) compared to Bølling-Allerød. A cooling YD subsurface is also simulated in the transient experiment of the global climate of the last 21,000 years (Liu et al., 2009(Liu et al., , 2012 Fig. 1) by applying meltwater south of 50 • N through St. Lawrence River (Fig. 1b), while H1 subsurface is warmed in response to meltwater discharge in between 50 • N and 70 • N NNA (Fig. 1a), consistent with many proxy records. Thus, the NNA subsurface temperature response seems to remain uncertain during YD and therefore requires special attentions.
The subsurface response in NNA has been suggested important in determining the oceanic feedback on the ice sheet calving (Shaffer et al., 2004;Bassis et al., 2017;Alvarez-Solas et al., 2010) at the onset of cold stadials as well as the recovery rate of the AMOC at the termination (Gong et al., 2013;Mignot et al., 2007;Liu et al., 2009). However, the mechanism for the subsurface temperature responses has remained unclear, because, conceptually, one may expect that a reduced northward ocean heat transport by the AMOC should induce a cooling in ocean subsurface as over the surface. Previous numerical experiments seem to suggest that the response of NNA subsurface temperature depends on the region of hosing: a hosing over the deep convection region of NNA generates a subsurface warming but a hosing south of NNA region favors subsurface cooling (Mignot et al., 2007;Brady and Otto-Bliesner, 2011), consistent with TRACE21 experiments (Fig. 1). The subsurface warming, under freshwater forcing over NNA, has been suggested to be caused by the suppression of deep convection, which prevents the subsurface water from mixing with the cold surface water (Mignot et al., 2007;Liu et al., 2009). In comparison, the mechanism for the subsurface cooling in response to freshwater forcing south of NNA is less clear. Some studies suggest that the cooling is caused by the active ventilation at the intermediate depth (Mignot et al., 2007;Liu et al., 2009;Ezat et al., 2014) while some other studies suggest that the cooling is caused by still active deep mixing (Brady and Otto-Bliesner, 2011). In all the cases, however, the freshwater forcing eventually reduces the AMOC and the associated northward oceanic heat transport; in the meantime, the freshwater forcing also reduces the mixing and deep convection. It therefore still remains unclear why the subsurface response depends critically on the region of hosing.
Here, we study the NNA subsurface temperature response by performing hosing experiments systematically in a fully coupled general circulation model. The mechanism of the subsurface temperature response is then investigated with a quantitative analysis of the temperature budget. The novel point of our study is to show that the subsurface temperature response in NNA is determined not by the location of the freshwater forcing, but by the effective freshening over NNA. A strong freshening, regardless of its forcing location, suppresses deep convection strongly in NNA and then induces subsurface warming, while a weak freshening, again independent of its location, mainly reduces the advective heat transport by the AMOC and therefore favors cooling.

CESM1.3 and experiments
The model employed in present study is the fully coupled National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM) version 1.3. CESM 1.3 is a state-of-the-art fully coupled atmosphere-ocean-land-sea ice general circulation model, with an atmosphere model using a finite volume dynamical core at a horizontal resolution of 1.9 • × 2.5 • and 30 vertical levels, and a primitive equation ocean model with a nominal 1 • horizontal resolution and 60 vertical levels. The horizontal resolutions of the land and sea ice models are the same as in the atmosphere and ocean, respectively. The climatology and variability of our current climate are properly reproduced by this model, and more details about other model components can be found in Hurrell et al. (2013).
The control experiment is driven by the forcing at the preindustrial level and is simulated for 500 years, reaching a quasiequilibrium state, and the last 100 years' output is used in our present study. Afterwards, following previous studies Brady and Otto-Bliesner, 2011;Otto-Bliesner and Brady, 2010), two sets of hosing experiments are branched from the control simulation, with the freshwater flux of various magnitudes applied over the NNA (50-70 • N) or the Gulf of Mexico (GOM, 15-33 • N, 80-105 • W), respectively. The two sets of experiments could be used as analogs for H1 and the YD cooling events, and here are denoted as NA# and GM#, respectively, with "#" indicating the hosing strength ranging from 0.1 to 1.0 for NA and from 0.5 to 2.0 for GM. For example, NA0.5 refers to a 0.5 Sv (1 Sv = 10 6 m 3 /s) meltwater flux evenly applied over NNA. It should be noted that the two sets of experiments employ the same set of boundary and initial conditions and only differ in the location and magnitude of freshwater forcing. All hosing experiments are integrated for 60 years when the NNA subsurface temperature anomaly becomes clearly visible (Fig. 2).

Temperature budget analysis
The temperature budget analysis is based on the thermal equation that governs evolution of sea water temperature in the CESM, It depicts the potential temperature tendency ∂θ ∂t is controlled by the zonal temperature advection (−u ∂θ ∂x ), meridional temperature advection (−v ∂θ ∂ y ), vertical temperature advection (−w ∂θ ∂ z ), heat source Q as well as horizontal M H (θ) and vertical M V (θ) mixing. Q includes heat source from ocean surface and solar penetration in the ∼100 m upper ocean, and is usually negligible in the subsurface ocean. Integrating C p ρ ∂Q ∂z from bottom to ocean surface leads to the net heat flux over ocean (Yang et al., 2015), where C p and ρ are sea water specific heat capacity and density. The vertical mixing, also often referred to as the diapycnal diffusion, is parameterized using the K-profile parameterization (Danabasoglu et al., 2006;Large et al., 1994) while the horizonal tracer diffusion represents the along-isopycnal diffusion based on the Gent-McWilliams parameterization (Gent and Mcwilliams, 1990). In the northern North Atlantic, over regions with vigorous deep convection (mixing layer depth >500 m), the vertical mixing is usually much larger than the horizontal mixing (Fig. S1f).
The annual mean output from the model is inserted into the temperature equation, such that advection, mixing and temperature tendency terms are reconstructed offline. This method reconstructs the advection term almost precisely in comparison with the online budget output in the intermediate-depth and deep ocean . Since the two mixing terms are highly time dependent and there are no standard outputs from the CESM, we combine these two terms as M(θ ) and calculate it as the residual of the equation. It should be noted that this method cannot be applied to the near-surface ocean where solar penetration plays additional role in the temperature budget expressed in Q, which is not presented in the model output.
For the subsurface, where Q is insignificant, the temperature transport equation is integrated with time from t = 0 to t = τ , such that where δ( ) means the difference between sensitive experiments and the control, and ( ) represents the time average τ 0 ( )dt τ between period from t = 0 to t = τ . This equation illustrates that the temperature anomaly at t = τ is the accumulation of advection and mixing over time, equivalent to the product of τ and averaged advection and mixing within the period. Unless stated otherwise, the main context uses the terms in square brackets for discussion.

NNA temperature response
The AMOC is reduced in both NA and GM experiments ( Fig. 2a) with the reduction more than 50% except for the two weakest forcing cases NA0.1 and GM0.5. Relatively, the AMOC weakens more rapidly in the NA experiments than the GM experiments of the same magnitude of freshwater forcing (e.g. NA0.5 vs GM0.5). However, given sufficiently strong freshwater in the GM experiments, the AMOC appears to be able to respond similarly to a NA experiment of a weaker forcing (e.g. NA0.5 vs GM2.0; Figs. 3b and d; Fig. 2a). As in previous studies, all experiments produce a robust surface bipolar seesaw response (Fig. S1) and subsurface warming in the tropical and South Atlantic (Fig. 3).
In contrast, the subsurface temperature response in NNA appears to be more complicated. The zonal mean temperature response suggests that the intermediate depth can exhibit both warming and cooling for hosing at either NNA or GOM (Fig. 3). For the strong hosing cases of GM2.0 and NA0.5 experiments, in spite of the dramatic surface cooling, a major warming anomaly develops around 500 m near 45 • N with a secondary warming in Greenland-Iceland-Norwegian (GIN) sea (Figs. 3b and d; Fig. S1). In the weak hosing cases of GM0.5 and NA0.1, however, the subsurface is cooled just below the surface cooling in the Labrador Sea (Figs. 3a and c; Fig. S1). Therefore, our hosing experiments show that the sign of NNA subsurface temperature response depends not only on the region of freshwater forcing, as suggested in previous studies (Mignot et al., 2007;Brady and Otto-Bliesner, 2011), but also on the magnitude of the forcing (Fig. 3). In particular, our GM2.0 experiment shows that freshwater forcing in GOM can also cause subsurface warming in NNA, if the forcing is sufficiently strong.

Advection vs mixing
We now examine the mechanism for the NNA subsurface temperature response quantitatively by diagnosing the temperature budget based on the thermal equation, with the focus on two competing effects generated by a surface freshening. A surface freshening over deep convection region reduces the mixing of warm subsurface and cold surface water masses, producing an anomalous subsurface warming (and surface cooling), while a surface freshening also weakens the AMOC, reduces the ocean heat transport into NNA and, in turn, cools the subsurface (and surface) water (Sgubin et al., 2017). Sign and magnitude of subsurface temperature response are thus determined by the competition between the two opposing effects, in contrast to the surface cooling coproduced by the reduction in ocean advection and mixing. We hypothesize that a subsurface warming is generated when convective mixing dominates advection while a subsurface cooling is generated when advection overwhelms convection. To verify this hypothesis, temperature budget analysis is performed at the subsurface NNA (40-70 • N, 300-1000 m). The region is selected based on the temperature responses in the hosing experiments.
It turns out that the climatology is maintained by the balance between the warming by total temperature advection of predominantly warm subtropical water (∼2 K/decade) and the cooling by local convective mixing of surface cold water (∼ −2 K/decade) (Table S1). With a strong freshwater flux over NNA (NA0.25 and NA0.5), the negative buoyancy flux suppresses deep convection rapidly in 5-10 years (as manifested in the March mixed layer depth) and causes subsurface warming (Figs. 4b and c) due to the reduced mixing (Table S1). The reduced mixing also traps cold water in the shallower surface mixed layer. Meanwhile, the weakening AMOC tends to cool subsurface owing to the reduced heat advection (Table S1). The strong warming due to reduced mixing overwhelms the cooling due to reduced AMOC advection (Table S1), leading to subsurface warming. This rapid suppression of deep convection (as seen in the rapid shallowing of mixed layer) and subsurface warming also occurs when a strong freshwater is discharged in GOM far away from regions of deep convection, as long as the forcing is sufficiently strong. Notably, this subsurface warming occurs in GM2.0, although this case exhibits an initial warming over the entire upper 2000-m in the first two decades (Fig. 4e). This initial warming is probably a result of the overshot of the AMOC and the associated advective warming, immediately after the freshwater discharge (Fig. 2a). In contrast to these strong hosing experiments, a subsurface cooling is generated when the hosing is weak, in either NA (NA0.1) or GM (GM0.5). The mixed layer shallows gradually over 50 years, giving rise to a subsurface cooling dominated by the advective cooling induced by the reduced AMOC.
The evolution of the subsurface temperature and the relative contribution of advection and vertical convective mixing can be shown more quantitatively in Fig. 5. The anomalous mixing is plotted against the anomalous temperature advection for each experi-ment (except for experiment NA0.1, whose response is too small to be visible; see Table S1 instead). All experiments start from the upper left corner and evolve downward, reflecting the advective cooling due to the reduced AMOC, and rightward, reflecting the warming associated with the reduced convection. Since the anomalous advection cooling and mixing warming are roughly balanced, each experiment evolves roughly along the diagonal line where the two terms are balanced exactly. The warming and cooling cases are found to be separated clearly by the diagonal line, with the warming cases above and cooling cases below, implying the dominance of mixing warming and advective cooling, respectively. This figure confirms our hypothesis and, moreover, quantitatively demonstrates that subsurface temperature response is determined precisely by the relative magnitude of advection and mixing, with a larger mixing for warming and larger advection for cooling. In comparison, hosing region is not the deterministic factor. A NA hosing, even of modest strength, usually leads to a subsurface ocean warming. If, however, the hosing and, in turn, mixing response, is too weak, subsurface temperature exhibits a cooling (such as NA0.1) (Table S1). A GM hosing usually leads to a cooling. If, however, the hosing and, in turn, mixing response, is sufficiently strong, it can lead to a subsurface warming (such as GM2.0). While the subsurface cooling generated by weak NNA hosing has been shown to be robust in previous works (Came et al., 2007;Stouffer et al., 2006), the large-scale subsurface warming generated by a strong GM hosing is, to our knowledge, to be the first time demonstrated in our GM2.0 experiment.

The effective freshening
It has been shown that, given a hosing region, NNA subsurface response tends to change from cooling to warming as the freshwater forcing is strengthened. The threshold freshwater flux from cooling to warming, however, differs for different hosing regions. In our model, for NNA hosing, the threshold is ∼0.2 Sv, as judged from the transition from cooling in NA0.1 to warming in NA0.25, while for GM hosing, the threshold is ∼1.5 Sv, as judged from the transition from cooling in GM1.0 to warming in GM2.0. It should be noted that the thresholds may change in longer sim-  x-axis is mixing anomaly (K/decade) and y-axis is total temperature advection anomaly (K/decade) including horizonal and vertical advection; shading corresponds to the temperature anomaly (K); the size of the dots corresponds to the time; the diagonal gray curve means advection is exactly balanced by mixing. See context for details. ulations as our simulations may not be equilibrated. The threshold freshwater flux is the minimum for hosing directly over the deep convection region of NNA, because freshwater over this region di-rectly suppresses convection before any dilution and is therefore the most effective. The meltwater discharged in remote regions, such as GOM, will be diluted before reaching the deep convection region in NNA and therefore requires a higher freshwater threshold for suppressing convection there. A natural question is therefore if there is an effective freshening forcing threshold that can be used to judge the temperature response, regardless of the hosing region.
One effective freshening threshold is the change of upper ocean stratification and, in turn, mixing intensity in the deep convection region over NNA. For hosing experiments, the change of upper ocean stratification is mainly caused by the surface freshening (Fig.  S2) due to freshwater forcing (NA cases) and transport (GM cases), with sea ice melting playing a modest role (Fig. S3) and vertical temperature change implying a destabilized effect (Fig. 3). The surface freshening in turn leads to anomalous salinity inversion over NNA (Fig. S2). In our model, once the salinity inversion exceeds a threshold of ∼ −4 psu, upper ocean stratification is stabilized, deep mixing is suppressed and the subsurface response changes from cooling to warming, regardless of the hosing region (Fig. 6a). Since mixing anomaly is largely linearly dependent on salinity gradient anomaly (Fig. 6b), the salinity threshold (−4 psu) corresponds to a mixing threshold of ∼1.2K/decade (Fig. 6b), which also appears to be the upper bound of the advection cooling (Fig. 6c).
The different responses seem to be also caused by a nonlinear relationship between mixing and AMOC reduction after hosing (δMix ∼ 1 β , β > 1) (see derivation in Supplementary material), which is likely caused by the nonlinear dependence of mixing on stratification and AMOC dynamics. This nonlinear mixing relation- ship competes against an almost linear relationship of advection and AMOC reduction (δ Adv ∼ ) (Fig. 6c). As such, advection overwhelms mixing to produce cooling for weak hosing, such as NA0.1, GM0.5 and GM1.0, while mixing overwhelms advection to produce warming for strong hosing, such as NA0.25, NA0.5, NA1.0 and GM2.0. The threshold for advective cooling is the same as that for mixing, of ∼1.2 K/decade (Fig. 6c). Thus, regardless of the hosing region, when the salinity inversion, and, in turn, mixing and advection changes over a threshold, NNA subsurface temperature response changes the sign.
While the nonlinear relationship between mixing and the AMOC is essentially determined by model parameters associated with isopycnal parametrization, the opposite responses of deep convection at Nordic Sea to the meltwater forcing of various magnitudes may also offer some insights to the relationship. Figure  S4 shows the spatial response of March mixing layer depth at the end of the simulation. Overall, the mixing layer depth shallows across the North Atlantic in all hosing experiments, with a dramatic shallowing over Labrador Sea. In the Nordic Sea, however, the response differs in strong hosing cases (NA0.5 and GM2.0) and weak hosing cases (NA0.1 and GM0.5), with strong (weak) hosing cases significantly shallowing (deepening). The deepening of mixing layer depth is likely linked to sea ice extension in weak hosing cases (Fig. S3) as the southward extension of sea ice leads to a positive salinity anomaly at the boundary of sea ice in the control run. In contrast, the shallowing mixing layer probably results from the applied freshwater forcing (Fig. S3), dominating the net salinity effect. As such the AMOC is largely shut off in the strong hosing cases and keeps active in the weak hosing cases (Fig. 2).

NNA subsurface response in other models
While our hosing simulations only last for 60 years that is far from true equilibrium as the adjustment time scale of the AMOC is usually longer than a few centuries Liu, 2013, 2014;Thomas and Fedorov, 2019;, the dependence of subsurface temperature responses on the effective freshening seems also robust in other climate models with longer hosing durations. The warming NNA subsurface with a strong effective freshening is revealed by a set of studies (Liu et al., 2009;Zhu et al., 2017;Mignot et al., 2007;Brady and Otto-Bliesner, 2011;Zhang et al., 2017), and the cooling NNA subsurface with a less effective freshening is also revealed from previous multi-model ensemble simulations (Came et al., 2007;Stouffer et al., 2006), in which 0.1 Sv freshwater discharges at the NNA directly. It is worth mentioning that hosing experiments with 0.35 Sv freshwater perturbation at NNA and subtropical North Atlantic are performed in an coupled climate of intermediate complexity, yielding a NNA subsurface warming and cooling respectively (Mignot et al., 2007). The sensitivity to hosing duration is examined in that model and fits our studies here.
The threshold of the effective freshening forcing from cooling to warming, however, may differ in other models. For instance, given a freshwater perturbation of 0.2 Sv over the NNA, Flückiger et al. (2006) produce a cooling subsurface, while the subsurface becomes warm for a 0.4 Sv perturbation. Recent studies by Swingedouw et al. (2013) and Thomas and Fedorov (2019) suggest that the NNA subsurface turns to be warm in response to 0.1 Sv freshwater perturbation, though multi-model ensemble simulations show a robust cooling at the same strength (Stouffer et al., 2006). In addition, the AMOC seems likely to partially recover after several hundred years due to the instability induced by the subsurface warming (Thomas and Fedorov, 2019). The conflicting responses may result from the hosing schemes. In Swingedouw et al. (2013) and Thomas and Fedorov (2019), the freshwater is applied around the coast of the Greenland to emulate a global warming melting scenario. This actually differs significantly from a "Heinrich hosing" where freshwater is usually evenly distributed at the North Atlantic (Stouffer et al., 2006) such that the freshening at the Labrador sea is likely to be less than that in the "global warming hosing". Therefore, the NNA subsurface becomes cooling in Stouffer et al. (2006) rather than warming in Thomas and Fedorov (2019). Nevertheless, the threshold of the effective freshening mostly relies on the nonlinear relationship between mixing and the AMOC, as advection is linearly dependent on the AMOC. More work are particularly encouraged to investigate the mixing-AMOC relationship by varying hosing schemes.
Finally, the warming subsurface revealed in our strong hosing cases (i.e. NA0.5) is a common feature in other peer studies (Liu et al., 2009;Zhu et al., 2017;Mignot et al., 2007;Brady and Otto-Bliesner, 2011;Zhang et al., 2017). The warming is expected as once the NNA surface freshens sufficiently, the deep convection fully stops and leads to a ∼2 K/decade warming in our model (Fig. 6c), but the entirely collapsed AMOC leads to a ∼1.5 K/decade cooling inferred from the intercept in Fig. 6c. The net effect is thus a ∼0.5 K/decade warming as shown in our NA0.5 and NA1.0 cases.
Physically, the inter-gyre connection (Mignot et al., 2007;He et al., 2019) between the subtropical and subpolar gyre contributes ∼0.5 K/decade energy input to high latitude North Atlantic in spite of the AMOC. As long as the deep convection disappears, the NNA subsurface becomes warm based on the background energy transport via gyres.

Implication to YD subsurface temperature response
Proxies suggest a modest freshwater discharge into the ocean via either the St. Lawrence River (Carlson et al., 2007) and Mississippi River (Kennett and Shackleton, 1975) or the Mackenzie River (Keigwin et al., 2018;Murton et al., 2010;Condron and Winsor, 2012) during the YD. These river discharges suggested by previous freshwater flux estimation (Carlson et al., 2007;Clark et al., 1999) and sea level reconstructions (Lambeck et al., 2014) are around 0.1 Sv, not likely to be the extreme GM2.0 case. This forcing strength is equivalent to a less effective freshening forcing at NNA regardless of where the freshwater discharges the Atlantic. As a result, it may trigger a cooling or, at least, less warming in the subsurface in YD relative to H1, as demonstrated by our numerical experiments. The cooling YD subsurface is also supported by the transient TRACE-21ka simulation (Fig. 1) in the Community Climate System Model 3, and one deglacial simulation in Parallel Ocean Program 2 (Zhang, 2016), the ocean component of CESM1.3.

Implication to trigger of Heinrich events
Our previous understanding about the trigger of the cooling events has been dominated by the assumption that freshwater discharge initiates the AMOC slowdown, reduces the northward heat transport, and leads to the Northern Hemisphere climate cooling, until recent high-resolution proxies reveal that these massive discharge occurs much later than the climate cooling and could hardly being the trigger (Barker et al., 2015). Instead, the massive freshwater discharges may be the results of the cooling events, with subsurface ocean warming associated with the AMOC variations being the key. The ocean warming increases underwater melt along the calving face of the Northern Hemisphere ice sheets and triggers rapid margin retreat and iceberg discharge and the following dramatic cooling (Shaffer et al., 2004;Marcott et al., 2011;Bassis et al., 2017;Ng et al., 2018).
However, it is so far unclear why Heinrich events occurred during some of the cold portion of millennial-scale climate oscillations called Dansgaard-Oeschger cycles but not all of them. The effective freshening mainly measured by the surface salinity change at NNA deep convection regions could be used as a practical indicator to quantify whether a freshwater signal, from subglacial melting and surface melting to iceberg discharges, should induce large-scale subsurface warming or cooling, and thus the "Heinrich" events during the Dansgaard-Oeschger cycles. Inversely, due to the largely linear dependence upon the effective freshening, the intermediate-depth temperature reconstruction over NNA could be also used to infer or verify the upper-level salinity change, thus the intensity of freshwater discharge as well as potential sea level change during these "Heinrich" events.

Conclusion
To summarize, our study quantitatively shows that the response of NNA subsurface temperature during Heinrich events is a result of the competition of local mixing associated with deep convection and advection associated with the AMOC. The sign and magnitude of the subsurface temperature depend on the effective freshening over the deep convection region, rather than where the freshwater is discharged as speculated in previous studies. A weak effective freshening favors the advective effect and causes subsurface cooling, while a strong effective freshening favors the mixing effect and generates subsurface warming.
There are some long term experiments that seem to be consistent with our experiments here (e.g. TRACE-21ka and Zhang, 2016), which lead us to believe that our conclusion may not be too sensitive to the length of experiment. It should be pointed out that our simulations here are short (decades-hundreds of years) compared with realistic paleoclimate events (thousands of years). Nevertheless, extended hosing simulations with varying convection parameters (e.g. vertical diffusivity) are highly desirable to further our understanding of the subsurface temperature responses on the effective freshening as well as the nonlinear relationship between mixing and the AMOC. It has been suggested that the greenhouse gas forcing may shallow the mixing layer and shut off the AMOC through air-sea surface heat flux exchange Gregory et al., 2005), thus the nonlinear relationship between mixing and the AMOC under anthropogenic warming is also a way worthy to be explored.

Data availability
The data supporting the findings of this study has been archived at http://dx .doi .org /10 .17632 /dtvvt9pc9c .2.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Computing and data storage resources, including the Cheyenne supercomputer (10 .5065 /D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR.