Comparison of Climate Model Simulations of the Younger Dryas Cold Event

Results of five climate model simulation studies on the Younger Dryas cold event (YD) are compared with a focus on temperature and precipitation. Relative to the Bølling-Allerød interstadial (BA), the simulations show consistent annual cooling in Europe, Greenland, Alaska, North Africa and over the North Atlantic Ocean and Nordic Seas with maximum reduction of temperatures being simulated over the oceans, ranging from −25 ◦C to −6 ◦C. Warmer conditions were simulated in the interior of North America. In two experiments, the mid-to-high latitudes of the Southern Hemisphere were also warmer, associated with a strong bi-polar seesaw mechanism in response to a collapse of the Atlantic meridional overturning circulation (AMOC). The modelled YD-BA temperature response was in general agreement with proxy-based evidence. The simulations reveal reduced YD-BA precipitation (up to 150 mm yr−1) over all regions with major cooling, and over the northern equatorial region. South of the equator, modelled precipitation seemed to increase due to a southward shift of the InterTropical Convergence Zone (ITCZ). The largest uncertainty in the YD is the high-latitude response, where the models show diverging results. This disagreement is partly related to uncertainties in the freshwater forcing. Most model studies assume an AMOC shutdown, but this is incompatible with proxy evidence.


Introduction
During the last termination, the Younger Dryas cold event (YD) interrupted the overall warming trend with the strongest response in the circum-North Atlantic region [1][2][3]. The YD started around 12.9 ka with an abrupt termination of the Bølling-Allerød interstadial (BA) and lasted until~11.7 ka, when a sharp warming signalled the start of the Preboreal (PB, Figure 1a). The abrupt climate changes associated with the YD have been the subject of many studies, based on both proxy evidence and climate model experiments [3,4]. In this paper, I evaluate the similarities and differences between five key modelling studies.
Over the last few decades, the YD has been the target of many modelling studies, using different model setups and experimental designs [5][6][7][8][9][10][11][12][13][14]. Most early work in the 1980s and 1990s was based on snapshot experiments performed with stand-alone atmospheric general circulation models (AGCMs), such as the GISS (Goddard Institute for Space Studies) model [10] or the ECHAM (European Centre-HAMburg) model [11]. In these snapshot experiments, the sea surface temperatures (SSTs) and sea ice conditions in the North Atlantic Ocean were prescribed, based on assumptions on surface ocean cooling and sea ice expansion. Later studies reported transient simulations with full dynamical coupling of the atmosphere-ocean system in which the ocean circulation was free to interact with the atmosphere [12,13]. These simulations relied heavily on the assumed forcing-history, in particular the freshwater forcing representing the melting of the ice sheets in North America and Europe. Essentially, the freshwater forcing was used in these simulations to tune the models to reproduce a climate evolution that is close to what is known from proxy records [12,13]. Recently, a different approach was used in which the model was constrained with selected proxy-based reconstructions in the North Atlantic realm through data-assimilation (DA [14,15]). In the latter approach, proxy data from the North Atlantic realm were used because the YD signal is most clearly expressed in this area [1][2][3], thus providing a clear constraint for the model.
Each of these three approaches has its merits. ACGM simulations are important to study the magnitude and distribution of the atmospheric response given specific changes in surface conditions (SSTs, sea ice, ice sheets, vegetation) (e.g., [16][17][18][19]). The transient coupled atmosphere-ocean simulations, on the other hand, are able to quantify how fast different climate system components respond to forcings and relative to each other, and how the response anomalies evolve in space (e.g., [20][21][22]). Finally, simulations with DA provide insights into what combination of forcings provides a climate response that is best fitted to proxy-based reconstructions and reveal what atmospheric and oceanic circulation is associated with the corresponding response (e.g., [23][24][25][26]).
However, there are also drawbacks associated with each of these three modelling approaches. First, AGCM simulations do not provide information on temporal changes in climate or on dynamical atmosphere-ocean interactions and imply strong assumptions on surface conditions. For example, AGCM simulations of the YD assume a collapsed Atlantic Meridional Overturning Circulation (AMOC), while there is no evidence for this in palaeoceanographic proxies [9,10]. Likewise, the scenarios applied in transient atmosphere-ocean simulations are not necessarily realistic. For instance, freshwater forcing scenarios may not be consistent with sea level records [12,13]. Similarly, the DA technique relies heavily on the proxy-based reconstructions, but these proxy data contain unresolved biases that are highly uncertain [25,26].
By exploring the similarities and differences in the results from these three approaches, insight is provided into the robust features of the YD and into the characteristics of climate variability during this important period. This is especially relevant for the fourth phase of the palaeoclimate modelling intercomparison project (PMIP4, [27]), with multiple modelling groups performing several climate simulations focusing on the last deglaciation 21 to 9 ka BP, including the BA-to-YD change [28].
In this paper, I therefore compare model results from five studies that cover the three discussed approaches: Rind et al. [10], Renssen and Isarin [29], He [13], Menviel et al. [12] and Renssen et al. [14]. I focus on the seasonal and annual response in temperature and precipitation, and address the following questions: -What is the range in temperature and precipitation response in the five models? -What robust features of the YD climate can be distinguished? -What aspects of the YD climate represent the largest uncertainties?

Models
In this paper, I discuss results from four different global climate models, being the GISS and ECHAM4 AGCMs, the CCSM3 fully coupled AOGCM and the LOVECLIM model of intermediate complexity. Table 1 summarises the main features.   [12], (c) ALL experiment reported by He [13], and (d) COMBINED experiment by Renssen et al. [14]. The temperatures in (b-d) are for the North Atlantic Ocean and Europe in the area between 60° W-30° E, and 40-70° N. Indicated are the timings of the samples for BA (small grey bar) and the YD (small black bar) discussed in this paper.

GISS Experiments
Rind et al. [10] performed two equilibrium experiments, both with boundary conditions for 11 ka (orbital parameters, ice sheet configuration and land-sea distribution). The timing of the YD start  [12], (c) ALL experiment reported by He [13], and (d) COMBINED experiment by Renssen et al. [14]. The temperatures in (b-d) are for the North Atlantic Ocean and Europe in the area between 60 • W-30 • E, and 40-70 • N. Indicated are the timings of the samples for BA (small grey bar) and the YD (small black bar) discussed in this paper.

GISS Experiments
Rind et al. [10] performed two equilibrium experiments, both with boundary conditions for 11 ka (orbital parameters, ice sheet configuration and land-sea distribution). The timing of the YD start at 11 ka was derived from the chronostratigraphy of Mangerud et al. [35] that is based on 14 C years, placing the YD between 11 and 10 14 C kyrs BP. The first experiment (11 kW) included modern sea Quaternary 2020, 3, 29 5 of 21 surface conditions and represents the climate during the BA interstadial. In a second experiment (11 kC) representing the YD, cold surface conditions were prescribed in the North Atlantic North of 25 • N following the CLIMAP [36] reconstruction for the last glacial maximum (LGM). Both experiments were run for 5 years, but only the output of the last 3 years was used for the analysis.

ECHAM4 Experiments
Results of two equilibrium experiments are discussed here (expGI1e and expGS1), with boundary conditions appropriate for the BA and the YD [29]. These boundary conditions include SSTs and sea ice cover, orbital forcing, ice sheet configuration, greenhouse gas concentrations, vegetation and land-sea distribution. Regarding the prescribed sea surface conditions, cooling and winter sea-ice cover were applied in the Nordic Seas in the BA simulation. In the YD experiment, strong North Atlantic cooling was prescribed [37], in summer based on foraminiferal evidence [38] and in winter on a modelled shutdown of the AMOC [39]. These two experiments were run for 12 years each, of which the last 10 years were used to account for a spin-up of the model.

CCSM3 Experiment
The results of He [13] discussed in this paper are from his ALL transient CCSM3 simulation, performed within the TraCE-21 ka project. The ALL simulation was started in the LGM and was run forward in time until 1990 CE, using transient forcings, including updates of ice sheets for every 500 years. Freshwater was added to the surface ocean at several locations around the North Atlantic and in the Ross and Weddell Seas. This applied freshwater forcing led to an AMOC collapse during the YD around 12 ka. From the transient results (Figure 1c), I selected two 100-year time periods representing the BA (13.1-13.0 ka) and the YD (12.1-12.0 ka).

LOVECLIM Experiments
The setup of the DG NS experiment by Menviel et al. [12] is similar to the ALL simulation reported by He [11], also starting in the LGM and using transient forcings. Also in this experiment, freshwater forcing was applied to the North Atlantic and the Southern Ocean, leading to an AMOC collapse in the YD, this time between 12.8 and 12.4 ka. Similar to the CCSM3 experiment, I selected two 100-year periods ( Figure 1b) to cover the BA (13.1-13.0 ka) and the YD (12.3-12.2 ka).
The COMBINED experiment of Renssen et al. [14] started with an equilibrium experiment for the BA, using appropriate forcings for 13 ka, including ice sheets, orbital parameters, greenhouse gas levels and freshwater forcing. This BA climate was used as a starting point for a transient simulation of the YD climate (Figure 1d), with a combination of forcings applied. The YD simulation was constrained through data-assimilation by proxy-based summer temperatures for Europe and selected annual mean SSTs, primarily from the North Atlantic region [14,40]. The YD results discussed here are averages of 96 ensemble members. Further details on the data-assimilation technique can be found in Dubinkina et al. [41]. The forcings consisted of a 3-year freshwater pulse of 5 Sv (1 Sv = 1 × 10 6 m 3 s −1 ) at the MacKenzie river outlet, a continuous meltwater forcing 0.1 Sv representing the background melt of the Fennoscandian and Laurentide icesheets, and a -0.52 Wm −2 negative radiative forcing, representing the reduction in atmospheric methane and nitrous oxide, and enhanced mineral dust levels. The applied data assimilation forced the model to an anomalous atmospheric circulation state with more frequent northerly atmospheric flow over Northern Europe. The freshwater forcing resulted in AMOC weakening by 50% relative to the BA, but not a shutdown. The recovery to the PB was obtained by turning off the data assimilation and by removing the background melt forcing. From this experiment, also two 100-yr periods were selected for the BA and the YD (Figure 1d).
To enable efficient analysis of all of these model results that were produced at various spatial resolutions, I converted all results to a common 1 • latitude by 1 • longitude grid, using the bilinear interpolation tool provided by Schulzweida [42]. Using these re-gridded results, I computed ensemble means, and, in addition, I calculated areal means for the change in temperature and precipitation from the BA to the YD. I focus on the following 15 areas: Greenland (70-80 •

Temperature
All models consistently simulate an annual mean YD-BA cooling in Europe and over the Nordic Seas and the North Atlantic Ocean (Figures 2 and 3). For these three areas, the models are also consistent in simulating a stronger temperature decline in boreal winter (December-January-February, DJF) than in boreal summer (June-July-August, JJA). Over Europe and the North Atlantic, the magnitude of this cooling is fairly coherent among the models, ranging from −2 to −6 • C for the annual mean YD-BA anomaly ( Figure 3). However, the values for the Nordic Seas show also the largest spread of all areas (Figures 3b and 4), where the LOVECLIM1.1 simulation suggests −19 • C for the annual mean, compared to around −2 • C for CCSM3 and ECHAM4. An analysis of Figure 2 suggests that this large deviation between models is related to differences in the implied latitude of North Atlantic Deep Water (NADW) formation. All simulations include strong YD-BA surface cooling over the mid-latitude North Atlantic Ocean, but the latitude of the strongest cooling depends on where the NADW formation takes place before the AMOC perturbation, and this latitude varies considerably among the models. For instance, CCSM3 has the main site of NADW formation between 50 and 60 • N, as reflected by the marked cooling (between −5 and −10 • C) between these latitudes in Figure 2d. In this CCSM3 result, however, the cooling in the Nordic Seas is much less (less than −5 • C), in contrast to the result of the two LOVECLIM simulations that have the main NADW formation site and strongest cooling (−10 to −25 • C) more northerly between 70 and 75 • N in the Nordic Seas (Figure 2c,e).
In Greenland, Alaska, Siberia, East Asia, and North Africa, four out of five models indicate YD-BA cooling for all the three timeframes. The deviating model is in most cases ECHAM4, with a warming signal during JJA. In Central North America, four models suggest warmer conditions in DJF during the YD relative to the BA, but LOVECLIM1.2 suggests colder YD conditions (−1 • C).
In the remaining regions, all in the Southern Hemisphere (SH), the responses are also diverging (Figures 2 and 3), depending more strongly on the experimental setup. In the two transient modelling studies that assumed a full AMOC collapse (CCSM3 and LOVECLIM1.1), positive YD-BA anomalies were simulated in most SH regions (+2 to +3 • C), most clearly in Antarctica, the Southern Ocean and the South Atlantic Figures 2c,d and 3k,o. Together with the marked cooling in the North Atlantic, this SH warming clearly reflects the bi-polar seesaw mechanism being active in these models [43][44][45]. When the AMOC is fully active, it transports important amounts of heat from the South Atlantic, over the equator to the North Atlantic Ocean [46], but when the AMOC collapses, this heat is accumulating in the SH, producing anomalous warming [47]. In contrast, the LOVECLIM1.2 study suggests on average negative YD-BA temperature anomalies in most areas of the Southern Hemisphere (Figure 2e). These relatively cool YD conditions reflect the impact of the negative radiative forcing (−0.52 Wm −2 , [14]), which in most SH areas overwhelms the warming impact of the relatively weak bi-polar seesaw in the LOVECLIM1.2 study [15].  Table 1 provides background information on the used models and on the experimental design.  Table 1 provides background information on the used models and on the experimental design.      Table 1.

Precipitation
The simulations indicate that most regions with YD-BA cooling experienced also a distinct decrease in precipitation (by as much as 100 mm yr −1 ), including Europe, North Atlantic, Nordic Seas, Greenland and North Africa (Figures 5 and 6). In areas further away from coldest region (i.e., the North Atlantic Ocean), such as Alaska, Siberia and East Asia, the precipitation response was more subdued and less consistent among the different models (Figure 6c,d,h). In Central North America, 4 out of 5 models simulate a modest increase in precipitation in the YD relative to the BA (+30 mm yr −1 ensemble mean, Figure 6g), in line with the warmer conditions suggested by most models. In the SH, ECHAM4 simulates negative YD-BA anomalies over most areas (Figure 6l-o), which is consistent with relatively cold YD conditions (Figure 2b). In contrast, the models with an active bi-polar seesaw (i.e., CCSM3 and LOVECLIM1.1) simulate mostly positive YD-BA anomalies in the SH (Figure 6j-o), in agreement with increased temperatures (Figure 2c,d).
The response in the tropics is generally similar to the pattern discussed by Renssen et al. [15], with negative and positive YD-BA precipitation anomalies north and south of the equator, respectively. This pattern is clearest in the simulations with full atmosphere-ocean coupling (LOVECLIM1.1/1.2 and CCSM3, Figure 5c,d) and is related to a southward shift of the intertropical convergence zone (ITCZ). This effect is well-known from other simulation studies with a perturbed AMOC, showing that the response enhances with a strong bi-polar seesaw effect (e.g., [13,46]). The two AGCM simulations (GISS and ECHAM, Figure 5a,b), with the surface forcing restricted to the North Atlantic region, do not show this contrasting North-South response in the equatorial region.

Precipitation
The simulations indicate that most regions with YD-BA cooling experienced also a distinct decrease in precipitation (by as much as 100 mm yr −1 ), including Europe, North Atlantic, Nordic Seas, Greenland and North Africa (Figures 5 and 6). In areas further away from coldest region (i.e., the North Atlantic Ocean), such as Alaska, Siberia and East Asia, the precipitation response was more subdued and less consistent among the different models (Figure 6c,d,h). In Central North America, 4 out of 5 models simulate a modest increase in precipitation in the YD relative to the BA (+30 mm yr −1 ensemble mean, Figure 6g), in line with the warmer conditions suggested by most models. In the SH, ECHAM4 simulates negative YD-BA anomalies over most areas (Figure 6l-o), which is consistent with relatively cold YD conditions (Figure 2b). In contrast, the models with an active bi-polar seesaw (i.e., CCSM3 and LOVECLIM1.1) simulate mostly positive YD-BA anomalies in the SH (Figure 6j-o), in agreement with increased temperatures (Figure 2c,d).
The response in the tropics is generally similar to the pattern discussed by Renssen et al. [15], with negative and positive YD-BA precipitation anomalies north and south of the equator, respectively. This pattern is clearest in the simulations with full atmosphere-ocean coupling (LOVECLIM1.1/1.2 and CCSM3, Figure 5c,d) and is related to a southward shift of the intertropical convergence zone (ITCZ). This effect is well-known from other simulation studies with a perturbed AMOC, showing that the response enhances with a strong bi-polar seesaw effect (e.g., [13,46]). The two AGCM simulations (GISS and ECHAM, Figure 5a

Comparison with Proxy-Based Evidence
I compare here the modelled temperatures with the global dataset of Shakun and Carlson [4], who provide proxy-based estimates of the annual mean temperature anomaly between the YD and

Comparison with Proxy-Based Evidence
I compare here the modelled temperatures with the global dataset of Shakun and Carlson [4], who provide proxy-based estimates of the annual mean temperature anomaly between the YD and BA. They used a range of proxies, primarily foraminifera-based Mg/Ca ratios and U k 37 , and δ 18 O on ice cores, foraminifera and speleothems. Shakun and Carlson [4] estimate that the global mean annual temperature was 0.6 • C lower in the YD relative to the BA. Two of the considered models (LOVECLIM1.2 and ECHAM4) simulate an equivalent global mean cooling of −0.6 • C cooling. In the other models in this study, the YD-BA global cooling was less than −0.6 • C, with −0.4 • C for GISS, −0.3 • C for LOVECLIM1.1 and −0.2 • C for CCSM3. The latter two results reflect the impact of a relatively intense SH warming associated with a strong bi-polar seesaw mechanism (compare Figure 2c,d with Figure 2a,b,e). On a global level, these relatively warm conditions in the SH compensate for a considerable part of the cooling in the circum-North Atlantic region. For a more detailed comparison to the individual proxy-based estimates of Shakun and Carlson [4], I divide the data into four zones of 90 degrees of longitude each and compare the proxy-based values with the model-based zonal means ( Figure 7). −0.3 °C for LOVECLIM1.1 and −0.2 °C for CCSM3. The latter two results reflect the impact of a relatively intense SH warming associated with a strong bi-polar seesaw mechanism (compare Figure  2c,d with Figure 2a,b,e). On a global level, these relatively warm conditions in the SH compensate for a considerable part of the cooling in the circum-North Atlantic region. For a more detailed comparison to the individual proxy-based estimates of Shakun and Carlson [4], I divide the data into four zones of 90 degrees of longitude each and compare the proxy-based values with the model-based zonal means (Figure 7).
In the Atlantic Ocean zone (Figure 7a), the models are relatively consistent with a marked YD-BA cooling (−6 to −8 °C) between 45 and 75° N that is in reasonable agreement with the proxy based estimates. An exception is the cooling simulated by LOVECLIM1.1 that appears about twice as strong as the rest of the data. In LOVECLIM1.1, this strong cooling results from a complete shutdown of the AMOC and subsequent cover by sea ice in the Nordic Seas [12]. Again, proxy-based reconstructions of the strength of the ocean circulation suggest that such an AMOC shutdown is unlikely during the YD [48][49][50][51].
The other zones show a divergence of the models at the SH mid-to-high latitudes (Figure 7b-d), suggesting a higher uncertainty here. LOVECLIM1.1 and CCSM3 show clearly warm YD conditions at SH high latitudes while the signal is less clear in the other models. The few available proxy-based estimates of YD-BA temperature change seem to support this warm SH anomaly, at least between 40 and 45° S [52][53][54][55][56]. Unfortunately, the dataset of Shakun and Carlson [4] does not provide estimates from the Southern Ocean around 60° S where the maximum YD warming is simulated. However, to the south of this modelled warming peak, Antarctic ice core records give clear evidence for warmer conditions during the YD time relative to the BA [57].  In the Atlantic Ocean zone (Figure 7a), the models are relatively consistent with a marked YD-BA cooling (−6 to −8 • C) between 45 and 75 • N that is in reasonable agreement with the proxy based estimates. An exception is the cooling simulated by LOVECLIM1.1 that appears about twice as strong as the rest of the data. In LOVECLIM1.1, this strong cooling results from a complete shutdown of the AMOC and subsequent cover by sea ice in the Nordic Seas [12]. Again, proxy-based reconstructions of the strength of the ocean circulation suggest that such an AMOC shutdown is unlikely during the YD [48][49][50][51].
The other zones show a divergence of the models at the SH mid-to-high latitudes (Figure 7b-d), suggesting a higher uncertainty here. LOVECLIM1.1 and CCSM3 show clearly warm YD conditions at SH high latitudes while the signal is less clear in the other models. The few available proxy-based estimates of YD-BA temperature change seem to support this warm SH anomaly, at least between 40 and 45 • S [52][53][54][55][56]. Unfortunately, the dataset of Shakun and Carlson [4] does not provide estimates from the Southern Ocean around 60 • S where the maximum YD warming is simulated. However, to the south of this modelled warming peak, Antarctic ice core records give clear evidence for warmer conditions during the YD time relative to the BA [57].
A major uncertainty is the YD signal over the Arctic Ocean (Figure 7). The simulations performed with LOVECLIM suggest cooler annual mean temperatures here (by 2 to 5 • C, Figures 2c,e and 7a-d), but the other models do not indicate any major cooling of more than 2 • C, and the CCSM3 result even shows warmer conditions over most of the Arctic (Figure 2d). Except for the Greenland ice cores that show marked cooling, the dataset of Shakun and Carlson [4] does unfortunately not include reconstructions of YD-BA temperature change for the Arctic. Other proxy-based evidence from lake and peat records from Eastern Siberia, Southern Alaska suggest distinct YD cooling [58] in agreement with the model results ( Figure 2). However, the YD warming reconstructed in Central Alaska, most of Northeastern Siberia and possibly Northern Alaska [58] is not reproduced by the models. Relatively dry YD conditions in Svalbard, as based on the extent of cirque glaciers [59], matches the simulated negative YD-BA precipitation anomaly in the Nordic Seas region ( Figure 5).
Concerning the seasonality of the YD-BA signal, all simulations suggest that cooling around the North Atlantic is stronger in boreal winter than in boreal summer (Figure 3a-f) in agreement with proxy data [60,61]. In Europe, the modelled YD-BA summer temperature anomalies (Figures 2 and 3e) agree well with 2 to 3 • C YD-BA cooling in summer as suggested by Renssen and Isarin [29] based on palaeobotanical proxies, and by Heiri et al. [40] based on chironomids. However, the models do not reproduce the proxy-based winter YD-BA cooling of at least 10 • C suggested by Renssen and Isarin [29], although the LOVECLIM1.1 result comes close (Figures 2 and 3e).
The modelling results discussed here give no indication for relative mild YD summer conditions in Europe recently suggested by Schenk et al. [9] based on both model simulations and palaeobotanical proxies. They applied a high-resolution (0.9 • × 1.25 • lat-lon) version of the CESM1 model that included coupled modules for the atmosphere-land-sea ice system to perform equilibrium simulations of the BA and the YD, using prescribed surface ocean conditions from the transient CCSM3 simulations [13] discussed in this paper. Schenk et al. [9] simulated similar summer temperatures in the YD and the BA in NW and SW Europe and even warmer YD summer conditions in Eastern Europe. Schenk et al. [9] find support for these relatively warm YD summers in reconstructions based on palaeobotanical proxies, using the climate indicator species method. The reconstructed warm summer conditions are based on an alternative interpretation of these palaeobotanical proxies relative to earlier work [29,62]. Further research on proxies for summer temperatures is therefore required to shed light on the nature of the summer YD climate in Europe.
The modelled large-scale patterns of precipitation change agree largely with the proxy-based evidence for hydroclimatic YD-BA change reviewed recently by Renssen et al. [15]. Over continents directly downwind from the cold North Atlantic Ocean, most notably in Europe and North Africa, the models simulate a much drier climate consistent with proxy-based reconstructions (e.g., [63][64][65][66][67][68][69]). However, proxies also indicate drying in East Asia (e.g., [70][71][72][73]), which is not clearly reproduced by the simulations (e.g., Figure 5f). In central North America, 4 out of 5 models suggest an increase in annual precipitation (Figure 6g) that is in agreement with different types of proxy evidence [74][75][76][77]. This enhanced precipitation is caused by an anomalous southerly atmospheric flow, bringing moist air from the Gulf of Mexico northwards to the American continental interior [15].
In the tropics, the models indicate a decrease in precipitation north of the equator and an increase to the south. This is a well-known model response to a weakening of the AMOC and the associated Northern Hemisphere cooling, causing the ITCZ to shift southwards [47]. Proxies from the tropics support this mechanism, showing for instance a wetter climate in Indonesia [78,79], but drying in the Philippines [79]. This North-South contrast in the hydroclimatic response is also found in proxy evidence from South America (e.g., [80]) and Southern Africa (e.g., [81,82]). In contrast to the precipitation, a major equatorial YD temperature signal is absent in both the models and the proxy data.
The proxy-model comparison of the result for LOVECLIM1.1 illustrates the uncertainties surrounding the forcing of the YD climate. On the one hand, the AMOC shutdown in this simulation produces a winter YD-BA cooling in Europe that is matching best with proxies, and a bi-polar seesaw that is consistent with proxy-based evidence for warming in the Southern Ocean and Antarctica. On the other hand, the imposed AMOC shutdown disagrees with palaeoceanographic evidence, and the associated relatively small global YD-BA cooling of −0.3 • C is a clear underestimation of the proxy-based value of −0.6 • C. This mismatch could indicate that the applied YD forcing is not realistic. A possible solution is the alternative YD forcing suggested by Renssen et al. [14], consisting of a weakened AMOC (~50%) with continued overturning in the Atlantic Ocean, producing a weak bi-polar seesaw with some warming at SH mid-to-high latitudes. This oceanic forcing would be combined with a sustained, marked negative radiative forcing related to lowered atmospheric levels of methane and nitrous oxide during the YD and an increased atmospheric dust load, for which ample evidence has been found in ice cores [34,[83][84][85]. In the LOVECLIM1.2 simulation, this negative radiative forcing results in YD cooling in many regions in the SH mid-to-high latitudes (Figure 2e), which is supported by some palaeoceanographic evidence [86]. However, the LOVECLIM1.2 result that includes this combined forcing does not produce a satisfactory match with proxy-based evidence everywhere either, as especially the European winter temperatures were underestimated.
The ability of models to reproduce proxy-based climate anomalies is also related to model-specific sensitivity to forcings, making it important to further investigate forcing scenarios with a suite of models as is planned within PMIP4 [28]. In this paper, I compared models of different set-up and complexity, and also the experimental designs were very different. In a future model-intercomparison it is important to apply a more standardised experimental design, and also to involve high-resolution models that include full atmosphere-ocean coupling to explore the inferences of Schenk et al. [9] for mild YD summers in Europe.

Conclusions
The comparison of model simulations shows consistent annual YD-BA cooling in Europe, the North Atlantic Ocean, Nordic Seas, Greenland, Alaska and North Africa. In all cases, the area of maximum YD-BA cooling is over the North Atlantic Ocean or the Nordic Seas, ranging from −25 • C to −6 • C. Except for North Africa, the temperature response in the tropics is minor in all simulations. Most models indicate YD-BA warming of up to 2 • C in the interior of North America. Over the Arctic Ocean and the mid-and high latitudes of the Southern Hemisphere, the models show limited agreement. A positive YD-BA anomaly of up to 4 • C was simulated over the Southern Ocean and Antarctica in the two coupled atmosphere-ocean experiments with an AMOC shutdown and a strong bi-polar seesaw effect.
The modelled YD-BA temperature change is generally in agreement with proxy-based evidence for maximum cooling in the circum-North Atlantic region and a strongest temperature response in the winter season. In the simulations discussed here, no indication is found for mild YD summer conditions in Europe, as recently suggested by Schenk et al. [9].
The simulated YD-BA precipitation change is clearly linked to the temperature response. In general, the climate became drier where it was colder, and wetter where it was warmer. The areas with a strongest YD-BA precipitation reduction (by up to 150 mm yr −1 ) are the North Atlantic Ocean, the Nordic Seas, Western Europe and the northern equatorial region. Over the southern equatorial region, a marked increase in YD-BA precipitation is evident, especially over the oceans. This general precipitation pattern is supported by proxy data.
The coupled atmosphere-ocean simulations with AMOC shutdown agree best with several aspects of proxy-based YD-BA climate response, such as the magnitude of winter cooling in Europe and the warming in the high-latitude Southern Hemisphere. However, the assumed full AMOC collapse is not supported by proxy data and these simulations strongly underestimate the magnitude of global YD-BA global cooling indicated by proxies. This mismatch indicates uncertainty in the YD forcing. It is proposed that negative radiative forcing could have been more important during the YD-BA climate change than assumed until now. This negative radiative forcing would represent higher atmospheric dust levels and reduced concentrations of methane and nitrous oxide during the YD, as evidenced by ice cores.
Funding: This research received no external funding.