Impacts of the PMIP4 ice sheets on Northern Hemisphere climate during the last glacial period

This study comprehensively investigates the impacts on the mean state of the Last Glacial Maximum (LGM) climate, particularly atmospheric circulation over the Northern Hemisphere associated with the different Paleoclimate Modelling Intercomparison Project Phase 4 (PMIP4) ice sheets, ICE-6G_C, GLAC-1D, and PMIP3, using the coupled atmosphere–ocean–vegetation model HadCM3B-M2.1aD. The simulation with PMIP3 ice sheets is colder than either of the two PMIP4 ice sheets mainly because of the larger area of land ice impacting surface albedo. However, changes in the circulation impact sea ice cover resulting in the GLAC-1D simulation being almost as cold. Although the PMIP4 ice sheets also induce different responses in the atmospheric circulation, some common features are identified in all simulations, including strengthening and lateral expansion of the winter upper-level North Atlantic jet with a large southwest-northeast tilt and summertime North Pacific jet, a southward shift of the wintertime Icelandic Low and Azores High and the summertime Pacific High. Compared to terrestrial-ocean reconstructions, all the PMIP4 ice sheet experiments overestimate the LGM cooling and wet conditions. The simulation with the ICE-6G_C ice sheet provides the closest reproduction of LGM climate, while the simulation with the PMIP3 ice sheet shows the coldest LGM climate state. Our study shows that in order to "benchmark" the ability of climate models to realistically simulate the LGM climate, we need to have reliable boundary conditions to ensure that any model biases are caused by model limitations rather than uncertainty about the LGM boundary conditions.


Introduction
The Last Glacial Maximum (LGM, ca. 21,000 years ago, 21 ka) is an important period in Earth's climate history. It was low atmospheric trace gases when large ice sheets were on northern North America, northern Europe, and parts of Siberia, and several ice sheets and ice caps were in the Southern Hemisphere. The climate at the LGM was generally colder and drier than today (MARGO Project Members 2009;Mix et al. 2001), with several ice sheets ice close to or at their maximal extent during this period (Clark et al. 2009). The LGM has been chosen as a target period for numerical climate simulations such as the Paleoclimate Modelling Intercomparison Project (PMIP; Braconnot et al. 2012;Kageyama et al. 2018), and numerous reconstructions have extensively documented the LGM environment and climate (e.g., MARGO Project Members 2009; Bartlein et al. 2011;Prentice et al. 2011;Albani et al. 2016;Cleator et al. 2020;Tierney et al. 2020a). It provides opportunities to compare and validate the climate simulations (e.g., Harrison et al. 2014;Kageyama et al. 2021) and test and improve the representation of specific climate processes and feedbacks in the models (e.g., Hopcroft and Valdes 2015). Furthermore, although LGM background conditions were different from the present, the globally averaged combined radiative forcing from LGM boundary conditions was similar in magnitude (of the order of 3-6 W m −2 ) to that expected for the end of the twenty-first century (Braconnot et al. 2012). The LGM has also been examined to constrain climate sensitivity (e.g., Rohling et al. 2012;Tierney et al. 2020b), an essential metric for future climate projections (e.g., Izumi et al. 2013;Harrison et al. 2015;Zhu et al. 2021), and the palaeoclimate simulations provide an opportunity to quantify systematic biases that are likely to be present in future climate projections. However, in order to "benchmark" the ability of climate models to simulate the realistic LGM climate, we need to have reliable boundary conditions to ensure that any model biases are caused by model limitations rather than uncertainty about the LGM boundary conditions. Boundary conditions are the factors affecting climate that are not simulated explicitly by a model; instead, they need to be prescribed. The crucial differences in the boundary conditions of the present compared to the LGM are (1) the extent and height of continental ice sheets, (2) insolation, and (3) atmospheric greenhouse gas (GHG) concentrations (Kageyama et al. 2017). Variations in continental ice sheets also cause changes in coastlines and bathymetry (i.e., land-ocean geography due to its link to sea level change), and in surface topography due to the viscoelastic properties of Earth's crust. For example, the expansion of the continental ice sheets in Greenland, Antarctica, northern North America, northern Europe, Siberia, Patagonia, and parts of Oceania at the LGM resulted in a eustatic sea-level lowering of ca. 130 m below the present level (Yokoyama et al. 2000;Peltier and Fairbanks, 2006;Lambeck et al. 2014). It brought changes in land-ocean geography, in particular, in Beringia, the Sunda Shelf in Southeast Asia, and the Sahul Shelf in Australia. Compared to the pre-industrial state, LGM insolation over the Northern Hemisphere decreased by up to 10 W m −2 in high latitudes during the boreal summer and increased by up to 3 W m −2 in mid-latitudes during boreal winter owing to changes in the Earth's orbital configuration in particular obliquity and eccentricity. However, overall, the change in the average global annual LGM orbital forcing was very small because of a similar precession.
The change in orbital forcing and greenhouse gas concentrations are well known (e.g., Berger 1978;Bereiter et al. 2015;Loulergue et al. 2008;Schilt et al. 2010). Similarly, the continental ice-sheet extents are relatively well constrained by geomorphological and geochronological data from moraines and glacial and fluvial deposits (e.g., Dyke 2004;Hughes et al. 2016), including limited evidence for the positions of ice streams and domes (e.g., Margold et al. 2018). However, while relative sea-level (RSL) records and Global Positioning System (GPS) measurements provide some information on the loading history (i.e., the spatial distribution of total mass by measurements of vertical and horizontal motion of the crust) of an ice sheet (e.g., Lambeck et al. 2014), there is little direct evidence for the icesheet thickness. Therefore, the estimation of LGM ice-sheet thickness has been implemented by combining numerical modeling and available data (e.g., Abe-Ouchi et al. 2015;Lambeck et al. 2010;Peltier et al. 2015;Tarasov et al. 2012), but the specification of ice-sheet topography has been a significant source of uncertainty in defining boundary conditions for LGM experiments (Kageyama et al. 2017). The uncertainty in the ice-sheet boundary conditions results from the different approaches used for the ice-sheet reconstruction (Ullman et al. 2014;Ivanovic et al. 2016).
Northern Hemisphere continental-scale ice sheets drive and amplify global climate change on different timescales (e.g., Clark et al. 1999). Several modelling studies have shown that the ice-sheet topography is likely to be the main factor altering the atmospheric circulation in the Northern Hemisphere (e.g., Kutzbach and Wright 1985;Manabe and Broccoli 1985;COHMAP Members 1988;Cook and Held 1988;Laîné et al. 2009;Rivière et al. 2010;Pausata et al. 2011;Löfverström et al. 2014;Ullman et al. 2014;Wang et al. 2018). Compared to the preindustrial simulations, CMIP-PMIP LGM simulations show that the North Atlantic westerly jets are stronger and narrower and extends into the North Atlantic (Beghin et al. 2016;Wang et al. 2018;Kageyama et al. 2021). Moreover, Wang et al. (2018) also point out a poleward shift of the LGM upper-level North Atlantic jet at the Northern Hemispheric scale, but most westerly jets rarely move except for the far-eastern North Atlantic where the shift of the westerlies' axis from peak to another. Through its effect on the storm tracks, changes in the position and the strength of the North Atlantic and North Pacific jets exert a strong influence on precipitation in western Europe and western North America, respectively (e.g., Oster et al. 2015;Beghin et al. 2016;Löfverström and Liakka 2016;Lora et al. 2017;Löfverström 2020). In some sensitivity experiments, a higher Laurentide Ice Sheet leads to a stronger and more zonally organized Atlantic jet, larger amplitude planetary waves, and Arctic warming (e.g., Held et al. 2002;Li and Battisti 2008;Likka and Löfverström 2018;Kageyama et al. 2021). Moreover, the seasonal-mean stationary wave response to the albedo and topography of ice sheets in the LGM were investigated; ice-sheet topography dominantly influences the wintertime stationary wave changes, while ice-sheet albedo mainly works for summertime stationary wave changes (Roberts et al. 2019).
The PMIP4 protocol for the LGM experiment provides three options of boundary conditions for the continental ice sheets: the ICE-6G_C, GLAC-1D, and PMIP3 reconstructions, and their associated topography and coastlines (Kageyama et al. 2017). The ICE-6G_C (VM5a) reconstruction combines a recent glacial isostatic adjustment (GIA) model and a global ice history model for reproducing the changes in total mass and ice thickness of major ice sheets (i.e., the Laurentide, Fennoscandian, Barents, Greenland, and Antarctic ice sheets) from the LGM to the present (Argus et al. 2014;Peltier et al. 2015). The total mass is constrained by a eustatic sea level that is well approximated by coral-based records. In contrast, the local variations of ice thickness are constrained by radiocarbon-dated RSL histories through the Holocene and by GPS measurements of the vertical and horizontal motion of the crust as well as GRACE satellite gravity measurements. The GLAC-1D reconstruction is combined from different sources (Tarasov et al. 2012;Briggs et al. 2014;Tarasov and Peltier 2002;Ivanovic et al. 2016), assembled for Ivanovic et al. (2016); The North American and Eurasian components are from Bayesian calibrations of a glaciological model, the Greenland component is a handtuned glaciological model, and the Antarctic component is from a scored ensemble of glaciological model runs. The four components employ ensembles of dynamical ice-sheet models that have been constrained with relative sea-level data and other available data (e.g., ice-core temperature profiles, geologically inferred deglacial ice margin chronologies, and present-day velocities and ice configuration). Finally, these components have been combined under GIA post-processing for a near-gravitationally self-consistent solution. The PMIP3 ice-sheet reconstruction (Abe-Ouchi et al. 2015) was created from three earlier individual ice sheet reconstructions for the LGM time slice (i.e., ANU, ICE-6G v2.0, and GLAC-1a) and allows direct comparison with the PMIP3 climate simulations that used these ice sheets (Braconnot et al. 2012). Refer to the literature cited above for detailed information.
This paper aims to investigate the impacts of each PMIP4 LGM continental ice-sheet topography on the Northern Hemisphere atmospheric circulation and resulting climate. It is the first study using all the PMIP4 ice sheets and focuses primarily on the large-scale robust changes as reliable climate responses to the ice sheet and inter-model differences resulting from slight variations in ice sheets. Although other studies have investigated varying LGM ice sheets, our study is the only one to quantify the impact of the CMIP6-PMIP4 ice sheets. Different modelling groups could use any one of these ice sheets benchmark their climate models, but there has been no systematic study of the uncertainties associated with these different ice-sheet reconstructions. Here we use one climate model, but the three different ice sheet models to understand and evaluate the range in the resulting climate. This paper proceeds as follows: in Sect. 2, we describe our model configurations and the experimental design. In Sect. 3, we evaluate the imposed forcing resulting from the change in ice-sheet and land-sea geography and surface albedo feedback for the LGM experiments. In Sect. 4, we assess the surface temperature changes and their decomposition in some target regions using a surface energy balance model. In Sect. 5, we describe how seasonal stationary waves, westerly jet, and precipitation patterns change in response to the LGM ice-sheet configurations. In Sect. 6, we evaluate the model performance with the PMIP4 LGM ice-sheet configurations through a palaeo data-model comparison. Finally, we conclude in Sect. 7.

Models
We use the coupled atmosphere and ocean model HadCM3 (Gordon et al. 2000), in particular, the version HadCM3B-M2.1aD reported by Valdes et al. (2017), which includes a relatively simple thermodynamic and free-drift sea ice model. The model configurations are similar to those in the original references but contain several small bug fixes and model improvements. The model also comprises an interactive dynamic vegetation model (TRIFFID) and the land surface scheme (MOSES2.1). The atmospheric component is run with 96 × 73 latitude-longitude grid boxes (3.75° × 2.5°) and 19 levels in the vertical, and the ocean component is run with 288 × 144 latitude-longitude (1.25° × 1.25°) and 20 vertical levels.

Experiments
For the control (0 ka) and LGM (21 ka) simulations, we have fully followed the experimental protocols for the CMIP6 piControl experiment (Eyring et al. 2016) and the PMIP4 LGM experiment (Kageyama et al. 2017). The LGM simulation includes atmospheric trace gases (CO 2 at 190 parts per million (ppm), CH 4 at 375 parts per billion (ppb), N 2 O at 200 ppb, with CFCs at 0), orbital parameters (eccentricity of 0.018994, obliquity of 22.949°, and perihelion-180° of 114.42°), and a solar constant that is the same as for the preindustrial (1361.0 W m −2 ). We derive our surface boundary conditions from the ICE-6G_C, GLAC-1D, and PMIP3 continental ice sheet reconstructions, and their associated topography, bathymetry, and coastlines.
LGM topography was first calculated as an anomaly from present-day topography on the original icesheet reconstruction grid and added to the HadCM3B's present topography after regridding onto the model grid.
We have followed the PMIP4 protocol for each of the ice sheets in order to ensure that we are testing the full differences between the choice of ice sheets. Hence the land-sea masks, bathymetry, and far-field surface elevations all change, as well as the ice sheet mask and ice sheet elevation. The PMIP4 protocol was developed to ensure that the imposed ice-sheet boundary conditions are geophysically consistent with associated changes in sea level and glacial isostatic adjustment. At the scale of the GCM, there are relatively small differences between the ice sheet boundary conditions for the land-sea mask (the difference between ice-sheet land-sea mask is less than 1% of the change between LGM and pre-industrial) and bathymetry (differences occur over less than 0.1% of the ocean surface and are never more than 1 model level). However, there are widespread but small differences in elevation over non-glaciated land (Fig. 1). This is partly due to the differing ice volumes resulting in different sealevel changes and partly due to the differing ice elevations (and history) resulting in differing crustal depression and rebound. We performed additional simulations in which the elevation changes beyond the ice sheet were unchanged from pre-industrial conditions. The resulting climates showed small temperature differences over land (consistent with lapse-rate estimated changes due to the ~ 100 m elevation change) but otherwise showed no statistically different responses. Hence our results do not depend on the elevation changes beyond the ice sheets.
Our gridded topography on the HadCM3B-M2.1aD model grid adequately captures the original LGM ice-sheet reconstructions. The North American ice sheet is composed of the Cordillera, Laurentide, Innuitian, and Greenland Ice Sheets, and the Eurasian ice sheet is composed of the Barents-Kara, Fennoscandia, and British-Irish Ice Sheets. However, there are substantial differences in the ice sheet configurations arising from each reconstruction, specifically regarding ice surface elevation and ice-covered areas (Fig. 1).
Compared to the GLAC-1D ice sheet, the ICE-6G_C and PMIP3 ice sheets has a taller Laurentide ice sheet, particularly the Keewatin Dome in west-central Canada and the Labrador-Quebec domes, and the elevation differences between the two ice-sheet reconstructions in the Northern Hemisphere are relatively small except the Cordilleran ice sheet (Fig. 1). In ICE-6G_C, the Keewatin dome extends towards the Cordilleran Ice Sheet and the relatively high-altitude Foxe dome. Compared to the other ice sheets, GLAC-1D provides the lowest ice-sheet surface elevation (less than 2500 m in the peak areas) across the entire North American ice sheet. In the European Ice Sheets, the ICE-6G_C and PMIP3 feature ice domes with around 2000 m on the Barents-Kara and Fennoscandia Ice . The values are differences between 21 and 0 ka using bright colours (ice grids) and pale colours (land grids). b the ice-sheet dif-ferences of ICE-6G_C and GLAC-1D from PMIP3. The ice-sheet mask is outlined in magenta Sheets, while the GLAC-1D reconstruction has a much taller dome on the Fennoscandia Ice Sheet with the maximum ice-sheet thickness above 2500 m. The Greenland ice sheet is similar among the three reconstructions. It is flatter (lower in the interior and higher at the margins) compared to the present-day configuration.
Both ICE-6G_C and GLAC-1D are the latest ice-sheet reconstructions for PMIP4 and have similar ice-covered areas over the Northern Hemisphere, while the PMIP3 ice sheets exhibit a much larger ice extent than the others, in particular over North America (Table S1). These differences can be important for the resulting climate since differences in the continental ice-sheet extent have significant impacts on surface albedo.
All simulations employed here were run for more than 500 years, initialised with a cold ocean from previous LGM simulations, and analyses are on the last 100 years of each simulation. For equilibrium climate states represented by the simulations, outgoing longwave radiation approximately balances the incoming absorbed solar radiation (i.e., the TOA net radiation values are close to zero; Table S2).

Shortwave radiative forcing and feedback analysis
To quantify the imposed forcing resulting from the change in ice-sheet and land-ocean geography and their associated feedbacks at the LGM, we used a simplified partial derivative method developed by Taylor et al. (2007). This method provides an accurate estimate of the large-scale radiative perturbation using standard model output. We build a simplified atmospheric shortwave radiation model to diagnose the planetary albedo (A) and its related three parameters [i.e., the surface albedo (α), the atmospheric absorption (μ), and the atmospheric scattering (γ)] from the simulated fluxes at the top-of-the-atmosphere (TOA) and at the surface (see Fig. 1 Taylor et al. 2007). These parameters differ between the LGM (21 ka) and control (0 ka) simulations. Changes in the net shortwave downward radiation (Q) at the TOA can be written as: where S represents the insolation, and Δ represents the difference between two climatic periods, 21 ka and 0 ka. The first term on the right-hand side represents the solar forcing induced by changes in the Earth's orbit. It corresponds to the difference in net solor radiation at the top of the atmosphere between 21 and 0 ka, if planetary albedo is the same. The second term results from changes in the planetary albedo, which is related to changes in α, μ, and γ. The planetary albedo changes can thus be decomposed in the following form: where e represents the nonlinear contribution of the different parameters and is negligible here. The effect of the change of each of the three parameters (α, μ, and γ) on shortwave radiation at the TOA can be thus estimated by computing the corresponding change in planetary albedo, keeping all but one parameter fixed to the values of the control simulation in Eq. (2).
In this study, changes in planetary albedo induced by changes in boundary conditions are considered as forcing. Thus, we defined surface albedo forcing as radiative responses from expanded LGM ice sheets and expanded land resulting from the lower sea level. On the other hand, changes to the surface or atmospheric properties resulting from the adjustment of the climate model to the new boundary conditions are considered feedback. Thus, we defined surface albedo feedback as radiative responses resulting from variations of snow cover, sea-ice cover, and vegetation.
The annual global averages of surface albedo forcing are the range of − 4.4 to − 3.4 W/m 2 (ice-sheet induced surface albedo forcing are − 4.0 to − 3.0 W/m 2 ) for the LGM ice sheets, the annual global averages of surface albedo feedback are − 3.7 to − 2.7 W/m 2 , and the global annual mean surface temperature changes are − 7.0 to − 5.6 °C. The PMIP3 ice-sheet experiment shows the strongest magnitude radiative responses and the coldest temperature (Table 1). The results are consistent with LGM simulations using other PMIP4 models (Kageyama et al. 2021). Compared to the previous studies with the PMIP3 ice sheet (e.g., Abe-Ouchi et al. 2015;Braconnot and Kageyama 2015), HadCM3B-M2.1aD with the PMIP3 ice sheet simulates similar surface albedo forcing (− 5.2 to − 3.6 W/m 2 ), but the temperature is 1 3 significantly colder (− 5.0 to − 4.4 °C in other models with the PMIP3 ice sheet ). Table 1 and Figure S1 shows seasonal (December-January-February, DJF; June-July-August, JJA) mean surface albedo forcings and surface albedo feedbacks over the Northern Hemisphere for the three LGM experiments. The DJF surface albedo forcings are much smaller in amplitude than those in JJA, especially in North America and northern Europe. This is because these regions are already mostly covered in snow at 0 ka, so the presence of an ice sheet under the snowy surface has limited effects on the surface albedo, and DJF insolation in high latitude is also very small compared to the JJA insolation. The area-weighted mean surface albedo-induced radiative responses (i.e., a sum of surface albedo forcings and surface albedo feedbacks) over the Northern Hemisphere range from − 4.5 to − 3.3 W/m 2 in DJF and from − 17.7 to − 14.2 W/m 2 in JJA. The smallest albedo radiative responses are simulated with ICE-6G_C and largest with PMIP3 ice sheets. These inter-model differences in surface radiative responses are caused by differences in the extent of the ice sheet as well as their effects on sea ice, and snow cover.
The surface albedo-induced radiative responses in boreal winter are small (− 10 to 0 W/m 2 ) over the Arctic but are more negative (− 70 to − 30 W/m 2 ) along the southern margins of the ice sheets ( Figure S1a). The DJF area-weighted mean surface albedo forcings over the Northern Hemisphere range from − 1.8 to − 1.3 W/m 2 among the LGM experiments, while there are large differences in surface albedo feedbacks (− 2.8 to − 1.9 W/m 2 ) mostly because of the different simulated sea-ice extents over the North Pacific and North Atlantic Ocean. Compared to the DJF sea-ice extent produced with ICE-6G_C, wintertime sea-ice cover in the GLAC-1D simulation is 235% and 145% greater in the Atlantic and the Pacific sectors (up to 60°N), respectively, and the sea-ice cover is 420% and 144% greater in the PMIP3 simulation.
The summertime area-weighted mean surface albedo feedback over the Northern Hemisphere is small and largely comparable (− 2.9 to − 2.3 W/m 2 ) between the LGM experiments, while differences in the surface albedo forcing are larger (− 14.8 to − 11.8 W/m 2 ) because of the continental ice-sheet extent variations ( Figure S1b). Compared to ICE-6G_C, the PMIP3 ice-sheet reconstruction has a 127% larger North American ice sheet and a 110% larger Eurasian ice sheet (Table S1). On the other hand, there are no large differences in the ice-sheet extent and their surface albedo forcing over the Northern Hemisphere between the ICE-6G_C and GLAC-1D. Over the continental ice sheets, the negative radiative forcing exceeds 100 W/m 2 in most places, and the larger negative values (less than − 150 W/m 2 ) are at the southern parts of the ice sheets. These features are consistent among the LGM experiments. Compared to the JJA sea-ice extent produced using ICE-6G_C, JJA sea-ice cover in the GLAC-1D simulation is 110% and 181% larger in the Atlantic and the Pacific Ocean, respectively, and with PMIP3 ice sheets, the sea-ice cover is 208% and 173% larger in the Atlantic and the Pacific sectors (up to 60°N). Figure 2a shows seasonal mean (DJF and JJA) surface temperature differences (21-0 ka) over the Northern Hemisphere. Area-weighted averages of surface temperature differences range from − 10.0 to − 8.0 °C in DJF and from − 7.0 to − 5.7 °C in JJA. Extreme wintertime cooling is simulated over the continental ice sheets in North America and Europe, over lands adjacent to the southern edge of the ice sheets, and oceanic regions typically covered by sea ice. Strong summer cooling is also simulated over the North American ice sheet, the south and central European ice sheets, and mid-latitude oceans, while summertime warming occurs over Beringia. This exceptional LGM warming in Alaska and East Siberia has also seen in proxy data and other paleo-simulations (e.g., Bartlein et al. 2011;Löfverström and Liakka 2016;Liakka and Löfverström 2018;Bakker et al. 2020;Cleator et al. 2020). The simulated warming happens in areas where ocean/sea ice becomes land because of the different heat capacity of land/ocean, including an increased seasonality in the regions, and changes in temperature advection and radiative properties. In order to diagnose the simulated temperature anomaly patterns, we use the same simple surface energy balance equation (detailed in the Supplementary Information) as one in Izumi et al. (2015), which was based on Laîné et al. (2009) and Lu and Cai (2009). The model deals with the radiative components (i.e., surface longwave and shortwave radiation) and non-radiative components such as surface sensible and latent heat fluxes, and the flow of heat into or out of storage for land and/or ocean (for oceans, this term includes the release of transported heat). These two radiative terms were also decomposed into the following five terms: the surface albedo effect (SAE), surface shortwave cloud forcing (SWCRF), surface longwave cloud forcing (LWCRF), surface clear-sky shortwave downward radiation, and surface clear-sky longwave radiation. Each of these terms shows the partial temperature change (PTC) contribution due to individual components of the energy balance to the total temperature anomaly, and the sum of these contributions will be approximately equal to the total surface temperature change. A possible reason for the less-than-perfect match is variations of surface emissivity, but we do not discuss it here. Figure 2b shows area-weighted averages of the seasonalmean surface temperature differences and the components of their partial temperature changes (PTC) over three regions: the Arctic (ARC: all grids between 60°N and 90°N), the North Atlantic Ocean (ATL: ocean grids between 40°N and 60°N in the Atlantic sector), and the North Pacific Ocean (PAC: ocean grids between 40°N and 60°N in the Pacific sector). Several components respond differently depending on the surface type (e.g., vegetated area, snow, ice sheet, sea ice, or ocean). Many of the responses of the components are robust, meaning the same sign of responses among the three LGM experiments, but the magnitude does vary.

The partial temperature change (PTC) analysis
Over the ARC, there are robust boreal wintertime and summertime surface temperature changes and their PTC, and those extents are similar among the three experiments. However, there are different responses of the PTC between the seasons; in DJF, changes to the surface clear-sky longwave downward radiation and heat storage mainly contribute to the surface temperature reduction, but changes to sensible heat and latent heat fluxes slightly moderate the cooling in the region. In JJA, the SAE and changes to the surface clear-sky longwave downward radiation decrease the surface temperature, but changes to the SWCRF and heat storage moderate cooling over the region. The sum of the non-radiative fluxes, including heat storage, contributes to surface cooling (− 2.3 to − 1.9 °C) in DJF and surface warming (7.8 to 9.3 °C) in JJA.
There are considerable differences in the simulated wintertime surface temperature changes over the ATL (− 19.3 to − 4.8 °C) and the PAC (− 12.2 to − 7.2 °C) due to wintertime sea-ice extent among the experiments. Over the ATL, changes to the surface clear-sky longwave downward radiation and sensible heat fluxes mostly contribute to the different amplitude of cooling among the experiments, whereas changes to heat storage induce strong warming. Over the PAC, changes to the clear-sky longwave downward radiation, SAE, and heat storage work for cooling, while latent heat fluxes strongly work for warming. Moreover, there are also large differences in the JJA surface temperature changes (− 8.6 to − 3.6 °C) over the ATL. The changes to the clearsky longwave downward radiation and SAE mainly work for cooling, but changes to SWCRF moderates this cooling. The other components contribute less to the surface temperature change. The summertime responses over the PAC are similar to the ATL, but smaller temperature differences (− 6.0 to − 4.4 °C) among the experiments. The sum of the non-radiative fluxes including heat storage contributes to wintertime surface warming over the ATL (1.1 to 2.1 °C) and over the PAC (0.5 to 0.8 °C), while it shows relatively small summertime responses over the ATL (− 1.2 to 0 °C) and over the PAC (− 0.6 to 0.1 °C).
In summary, the surface clear-sky longwave downward radiation is one of the crucial components for the surface temperature changes over the target regions during both seasons. In JJA, the two other components also show robust responses: SAE drives cooling over the ARC, and a SWCRF change drives warming. In addition, changes to non-radiative fluxes including heat storage show substantial contributions to the surface temperature, in particular the DJF. The non-radiative fluxes largely contribute to summertime surface warming over the ARC but smaller responses over the ATL and PAC. However, their response varies depending on regions, seasons, and experiments.
As mentioned in Sect. 3, the HadCM3B-M2.1aD LGM simulation with the PMIP3 ice sheet has significantly lower surface temperatures than in other models employing the PMIP3 LGM ice sheets. Using this surface energy balance model, we further investigate the difference of surface temperatures among the PMIP3 LGM simulations. We find that compared to the PMIP3 ensemble mean data , our global mean (85°S to 85°N) surface temperature is more than 2 °C colder, mainly because of clear-sky longwave downward radiation changes during the winter ( Figure S2). This change in the surface clear-sky longwave downward radiation is a result of the feedback between the enhanced cooling and changes in atmospheric water vapour, atmospheric CO 2 , and atmospheric energy transport.

Atmospheric circulation
This section investigates the responses of Northern Hemisphere atmosphere circulation, especially with regards to stationary waves, low-and high-level westerlies, and precipitation, to the continental ice sheets at the LGM during boreal winter (DJF) and summer (JJA).

Evolution of the winter circulation and stationary waves
In the control state, the wintertime sea level pressure (SLP) field over the Northern Hemisphere can be broadly characterised as subpolar low-pressure systems and subtropical high-pressure systems over oceans (Fig. 3a). The Aleutian Low and the Icelandic Low are in the North Pacific and the North Atlantic, respectively, and the North Pacific High and the Azores High are to the south of each low system. The Siberian High lies over central Asia. Compared to the control simulation, the LGM simulations increase the winter SLP over the Northern Hemisphere ( Fig. 4a; Figure S2). At the LGM state, large-scale SLP responses are mostly influenced by topography (Pausata et al. 2009(Pausata et al. , 2011, and one of the robust features of the response is a southward shift of the climatological mean SLP centres relative to the PI simulations (e.g., Pausata et al. 2011). During the LGM winter, compared to the control state, the Icelandic Low and Azores High centres move south in all our simulations, and a southward shift of the Aleutian Low centre happens in the ICE-6G_C and PMIP3 simulations ( Fig. 3a; Table S3). Moreover, the Icelandic Low centre moves eastward, particularly in the PMIP3 simulation, and the Azores High centre robustly moves westward, particularly in the ICE-6G_C simulation. Compared to the control state, the LGM simulations also develop stronger winter SLP gradient in the North Atlantic and Pacific sectors, while changes in the winter SLP gradient between the Aleutian Low and Siberian High are inconsistent among our LGM experiments (Table S4). Figure 3b shows the climatological boreal winter 850 hPa and 200 hPa eddy geopotential height (an eddy is defined here as a deviation from the zonal mean) in the control state, revealing several troughs and ridges. HadCM3B-M2.1aD does not adequately simulate the low-and high-level deep troughs from northwest Atlantic to northeast Canada shown in ERA-interim data (Dee et al. 2011), which bring cold Arctic air into central and eastern North America (Fig. 1a in Roberts et al. 2019). However, the model simulates realistic wintertime low-level and upper-level stationary waves in both the large-scale patterns and their amplitudes.
There are a set of troughs and ridges over the mid-and high-latitudes in the control state; low-and high-level troughs over the western North Pacific extending eastern Asia and Siberia and low-and high-level ridges over western North America and the eastern North Atlantic extending over Europe (Fig. 3b). The horizontal trough-ridge patterns are displaced westward with increasing height (i.e., westward tilted) because of meridional temperature advection by the associated geostrophic wind, which induces cooling to the west of the low (Holton 2004). Vertically coherent structures of the ridges are over the eastern Atlantic associated with Azores High, while the low-level ridge over Eurasia is associated with the Siberian High, but there is no corresponding ridge at the upper level. Figure 3b also show the wintertime stationary waves for the LGM experiments (the LGM anomalies from the control state are shown in Figure S3). The large-scale spatial patterns are similar in each simulation, but there are notable differences in the wave patterns between the control and LGM experiments over the North America and downstream/upstream of the Laurentide ice sheet. The upper-level deep trough is over the central to the east part of the North American ice sheets (i.e., the Keewatin-Laurentide ice sheets) and extends poleward-downstream of the Laurentide ice sheet, especially in the GLAC-1D and PMIP3 experiments.
There are large differences between low-and high-level troughs over the northeast Atlantic, downstream of the Laurentide ice sheet; the ICE-6G_C experiment shows stronger low-and high-level ridges over Greenland and the Denmark Strait, while the GLAC-1D and PMIP3 experiments show low-and high-level deeper troughs. Another notable difference is an upper-level strong ridge over the equatorwarddownstream of the Laurentide ice sheet and over Europe to the Tibet plateau. On the other hand, the upper-level ridge in the LGM simulations shifts to Alaska and Bering Strait, upstream of the Laurentide ice sheet, and the lowand upper-level ridges are much stronger, especially in the PMIP3 experiment. All experiments also simulate low-and high-level troughs over the central-eastern North Pacific. However, compared to the control state, deeper troughs at the central area are shown in the ICE-6G_C experiment, while shallower troughs throughout the North Pacific midlatitudes are shown in the GLAC-1D experiment. The prevailing westerly wind is prominent in the middle-latitude troposphere. First, we focus on the zonal winds at the 200 hPa level, featuring the maximum speed of westerly winds in the control simulation. Figure 4a show the boreal winter 200-hPa westerly jets in the control and LGM climates. The upper-level westerly wind in the Northern Hemisphere is zonally asymmetric with discontinuous areas of the wind field over the eastern North Pacific and North Atlantic due to the land-ocean distribution. Compared to the control simulation, the LGM simulations have robust features (Fig. 4a, Figure S4); decreases in the zonal wind speed maxima at 200 hPa (around 5-20% weakening) and their locations to the north at the Northern Hemispheric scale. However, 200 hPa westerly jet over the North Atlantic alone strengthens and extends to central North Atlantic, and its core moves poleward. Because cores of jets move little in other areas, the 200 hPa wintertime North Atlantic jet in the LGM climate can be interpreted as tilting heavily to the southwest-northeast rather than moving poleward across the Atlantic jet (Fig. 4a). The strengthening of the North Atlantic jet is consistent with recent PMIP simulations (e.g., Löfverström 2020; Kageyama et al. 2021), while the jet with a large southwest-northeast tilt is consistent with the recent simulations (Löfverström 2020; Figure S5) except PMIP1 simulations (Li and Battisti 2008). Figure 4b shows the boreal winter 850-hPa westerly jets in the control and LGM climates. The 850 hPa atmosphere characterizes the eddy-driven jet and its fluctuations while remaining above the surface layer. Compared to the surface westerlies, the 850-hPa zonal wind exhibits accordant behaviour and excludes the impact of the boundary layer processes (Chavaillaz et al. 2013). The control simulation shows the strong westerlies over mid-latitudes of both the North Pacific and North Atlantic Oceans. Compared to the control state, all the LGM experiments show the stronger westerlies across the mid-latitudes of the North Atlantic, but the ICE-6G_C experiment has a different spatial pattern, the low-level jet with a zonalization, while the other simulations show jets with larger southwest-northeast tilt. The ICE-6G_C experiment also shows the strengthening and equatorward shift of the low-level Pacific jet.  Figure 4 also shows the wintertime LGM precipitation anomalies from the control state. The wintertime areaweighted average of precipitation over the Northern Hemisphere is 2.4 mm/day for the control state, and precipitation within the primary Atlantic and Pacific storm tracks, which display a southwest-northeast tilt, reaches more than 7 mm/ day ( Figure S4c). Precipitation along the Pacific coast of North America is largely controlled by North Pacific storm tracks, which are manifested in the position and intensity of the Aleutian Low (Salathé 2006;Oster et al. 2015). Wintertime precipitation over the western Mediterranean is mainly influenced by the North Atlantic westerlies associated with the jet stream position (Cortesi et al. 2014;Beghin et al. 2016). Compared to the control state, our LGM experiments show 13 ~ 17% decreases in area-averaged precipitation over the Northern Hemisphere, especially the storm track regions and the continental ice-sheet areas. The ICE-6G_C experiment shows the precipitation decrease due to shifting the precipitation area to the south of the Laurentide Ice Sheet with winter westerly storm track moving south (Fig. 4b).
The feature has also been reported in previous studies (e.g., Oster et al. 2015;Löfverström 2020;Tulenko et al. 2020). The two other LGM experiments show the precipitation decrease without a southward shift of the westerly storm track, possibly due to the thermodynamic mechanism. On the other hand, all the LGM simulations show an increase in precipitation over the Iberian Peninsula and the west coast of Europe. The GLAC-1D and PMIP3 experiments show more precipitation along the northern European coast because of the heavier tilted and stronger North Atlantic jets.

Evolution of the summer circulation and stationary waves
In the control state, the summertime SLP (Fig. 5a) is almost a reverse of the wintertime SLP: thermal lows on the continents and subtropical highs over oceans. The Bermuda (i.e., Azores High in winter) is over the Atlantic Ocean and the Pacific High over the Pacific Ocean, while the Icelandic Low and Aleutian Low retreat northward during summer. Compared to the control experiment, the LGM experiments show that the Bermuda High shifts northeast and Pacific High shifts southeast, respectively ( Fig. 5a; Table S3). While wintertime stationary waves extend and peak further poleward (Fig. 3b), summertime stationary waves are primarily found in the NH subtropics (Fig. 5b). Compared to the wintertime stationary waves (where the role of orography is large), in the summertime, diabatic heating, which is closely related to land-ocean thermal contrast, plays a much larger role in maintaining the subtropical stationary waves, because the zonal wind speed is weaker in summer than in winter (Ting 1994). The HadCM3B-M2.1aD model simulates spatially realistic summertime stationary waves, but their amplitudes are slightly stronger than reanalysis data suggests they should be (Roberts et al. 2019). The control simulation shows a series of baroclinic structures with surface ridges and upper-level troughs over the Pacific and Atlantic Oceans and with low-level troughs and upper-level ridges over central Asia and southwestern North America (Fig. 5b).
Compared to the control state, all the LGM experiments show several robust responses to large-scale summertime stationary wave: a stronger low-level ridge over North America, a weaker upper-level ridge at southwest North America, and a deeper upper-level trough at northeast North America ( Fig. 5b; Figure S3). These simulations also show a weaker low-level ridge and a shallower upper-level trough over North Pacific, and stronger low-and upper-level ridges at the North Atlantic and western to central Eurasia. Responses of summertime stationary waves at the LGM can be interpreted as the combined responses to the reduced diabatic heating and the mechanical forcing of the raised topography (Roberts et al. 2019). The change in summertime surface albedo relating to the prescribed ice sheets is much more important than wintertime, partly because a major part of northern North America, Scandinavia and Siberia is widely and thickly snow-covered in winter, and partly because the albedo has a much stronger effect during summer when the solar forcing is stronger. Hence, the summertime stationary wave responses are strongly influences by the ice-sheet area and more consistent than ones in winter.
Upstream of the North American Ice Sheet, the summertime stationary wave responses can also be interpreted as the response to the implied heating anomaly (Fig. 5b; Figure  S3). The summertime heating over North America strongly influences subtropical high over the North Pacific (Rodwell and Hoskins 2001). Reduced heating over North America weakens the low-level subtropical high and the upper-level trough over the Pacific, showing the strongest response in the PMIP3 ice-sheet experiment. On the other hand, downstream of the North American Ice Sheet, some ridges extend throughout the atmosphere with a minimal tilt over the North  Figure 6a shows the boreal summer 200 hPa westerly jets in the control and LGM states. Compared to the wintertime jets, the summertime zonally averaged mid-latitude jets are much weaker and shifts poleward. Compared to the control experiment, all the LGM experiments have several robust responses; North Atlantic and North Pacific jets rarely move latitudinally, but the North Pacific jet strengthens and extends across the mid-latitude Pacific ( Figure S6). Moreover, the boreal summer 850 hPa westerly jets in the control and LGM climate states are in Fig. 6b. Compared to the control state, all the LGM simulations show that the North Pacific jet strengthens with a larger southwest-northeast tilt and shifts to the equator. Finally, Fig. 6 also shows the boreal summer precipitation anomalies from the control state. Compared to the control experiment, all the LGM experiments have some robust responses: (1) more precipitation at the west-central USA, possibly due to strengthening of North Pacific jets and their equatorward shift and (2) less precipitation over eastern North American and some monsoon regions such as West Africa, South Asia, and East Asia.

Data description
To assess the realism of the simulated large-scale climate change over the Northern Hemisphere at the LGM with a given ice sheet, we compare the model output with existing palaeodata syntheses; pollen-and plant macrofossilbased climate reconstructions from Bartlein et al. (2011) and surface ocean reconstructions from the MARGO data set (MARGO Project Members 2009) and from Schmittner et al. (2011). These data sets provide a grid-cell mean of reconstructed anomalies (i.e., 21-0 ka) and a grid-cell mean of standard errors of estimates on a regular 2 × 2° grid (Figure S4 and S5). The grid-cell value was estimated by simple averaging using individual reconstructions within the grid, and the uncertainty was obtained as a pooled estimate of standard error calculated in the usual way. The gridded land data is the same as in Izumi and Bartlein (2016) and involves a simple area average of all reconstructions within a 125 km radius search window centered at each grid point. The gridded ocean data is the same as described in Harrison et al. (2014). The paleoclimate data sets include several sources of uncertainty (e.g., age modelling and the reliability of reconstruction approaches and calibrations). In creating these data syntheses, both sites with poor chronologies and records from atypical environments were removed to consider such uncertainties Harrison et al. 2014).
Moreover, we used bioclimatic variables from Cleator et al. (2020), which used 3D variational data assimilation techniques with a prior derived from the PMIP3-CMIP5 LGM simulations and terrestrial reconstructions from the Bartlein et al. (2011) and other data. Cleator et al. (2020) provides bioclimatic data on the 2 × 2° grid at more points than original terrestrial reconstructions and extends the geographic coverage (e.g., eastern Siberia, Figure S8 and S9).
To evaluate the simulated annual temperature and precipitation and the simulated seasonal temperatures, we selected seven available climate and bioclimatic variables: mean annual temperature (MAT), total annual precipitation (TAP), mean temperature of the coldest month (MTCO), and mean temperature of the warmest month (MTWA) over land grids; boreal summer (July-August-September, SSTsum), winter (January-February-March, SSTwin), annual (SSTann) sea surface temperatures over ocean grids. Regarding the target areas, we also selected seven land areas and two ocean areas with abundant reconstructions and bioclimatic data over the Northern Hemisphere: Western Europe

Methodology
To provide a primary measure of large-scale agreement between observations and model output, we use areaweighted averages by the area of 2° × 2° regular grid cells (Fig. 7). The area-averaged bioclimatic variable changes over each target region were calculated with only the grid cells that include paleoclimatic data. For the approach, the model output was interpolated to a regular 2° × 2° grid using bilinear interpolation. The uncertainties associated with the reconstructions (i.e., prediction error variances) and the use of different methodologies were propagated into the area averages using a Monte Carlo approach (Izumi et al. 2013) as detailed below.
For the Bartlein et al. (2011) and MARGO (2009) reconstructions, the error bars include reconstruction uncertainties and bootstrap estimates of the spatial means' standard error (Fig. 7). The reconstruction uncertainties were calculated as follows: for any one weighted average (whole or one of the bootstrap iterations), the reconstruction uncertainties were included by generating a random number from the normal distribution for each reconstructed value and then rescaling that value using the standard error of the reconstruction. We used this approach to generate 1000 replicate estimates of each reconstructed value. Thus, the bootstrap estimates of the standard error of the reconstruction include the variability due to the particular sample of points included in the mean, the reconstruction uncertainty of those points, and the overall spatial variability.
The error bars for each simulation and the Cleator et al. (2020) data in Fig. 7 are bootstrap estimates (1000 replications) of the spatial means' standard error the data locations.
For each replication, the values are sampled with replacement, and the area-weighted mean is calculated. Thus, the variability among the 1000 values is related to the particular sample of points selected at each iteration. The error bars (two standard deviations here) show the model output's uncertainty related to the observations' particular distribution. MAT mean annual temperature, MTCO mean temperature of the COldest month, MTWA mean temperature of the WArmest month, TAP total annual precipitation, SSTann annual mean SST, SSTwin winter mean SST, SSTsum summer mean SST

Results
There are no notable differences between the Bartlein et al. (2011) andCleator et al. (2020) data sets for most of the temperature variables across their common grid cells (Kageyama et al. 2021). However, there are large differences in area-weighted averages of some bioclimatic variables in some regions because of the different spatial homogeneity (i.e., low density and uneven distribution of the Bartlein et al. (2011) data): ΔMAT in western Europe (5.2 °C), eastern Europe (4.2 °C) and Alaska (3.7 °C) and ΔMTCO in eastern Europe (3.8 °C) and western North America (3.9 °C). ΔTAP in western North America (302 mm) and Alaska (163 mm) also show a large difference between the two data sets. The features of the differences are as follows: compared to the Bartlein et al. (2011) data, the Cleator et al. (2020 data show (1) much lower annual and winter temperatures in Europe, (2) much lower annual temperature and more precipitation decrease in Alaska, and (3) more precipitation in western North America.
Features of large-scale data-model comparison over lands are described here (Fig. 8). Compared to the Bartlein et al. (2011) data, the bioclimatic variables and regions for which at least one model results match the reconstruction are as follows: • ΔMAT in western Europe, eastern Europe, and eastern North America • ΔMTCO in eastern Europe, western North America, and eastern North America • ΔMTWA in Alaska and eastern North America • ΔTAP at western Siberia, western North America, and eastern North America.
Simulated bioclimatic variables by the ICE-6G_C experiment or both ICE-6G_C and GLAC-1D experiments are put into the reconstructions ranges at most cases above.
Compared to the Bartlein et al. (2011) data, all the simulations largely overestimate annual and seasonal cooling over North America, in particular in Alaska: annual cooling (− 12.0 to − 8.8 °C) and winter cooling (− 10.3 to − 5.1 °C). Annual cooling in West Siberia (− 6.5 to − 4.5 °C) and summer cooling western Europe (− 7.0 to − 3.9 °C) are also overestimated. The ICE-6G_C experiment shows the least cooling (i.e., closest to the reconstructions) at the target regions above, while the PMIP3 simulation shows the largest cooling at the target regions above except for Alaska, where the GLAC-1D simulation shows the coldest temperature. On the other hand, all simulations largely underestimate winter cooling in western Europe (10.2-11.1 °C), regardless of the differences of ice sheets and extension of sea ice. Compared to the reconstruction, all simulations also show smaller precipitation decreases in western Europe (425-448 mm) but much larger precipitation decreases in Alaska (− 341 to − 310 mm).
Compared to the MARGO data set, the ICE-6G_C simulation shows the appropriate SST decrease over the midlatitude North Atlantic (32°N-60°N), while the other LGM simulations overestimate the cooling with large extension of sea ice. All the LGM simulations slightly overestimate cooling over the high-latitude North Atlantic (60°N-80°N).

Discussion
Compared to the reconstruction, all the LGM experiments show an adequate change in the annual mean temperature changes in Europe, especially in the western region. However, those seasonal temperature changes are very poor: underestimating winter cooling and overestimating summer cooling. This underestimating winter cooling has long been pointed out in the LGM simulations with other climate models (e.g., Braconnot et al. 2012;Kageyama et al. 2021), while the change in simulated sea surface temperature over the mid-latitude North Atlantic does not show an underestimation of this cooling. On the other hand, all the LGM experiments show similar total annual precipitation to the control one (i.e., underestimation of precipitation changes). Ullman et al. (2014), Beghin et al. (2016), and Löfverström (2020) show increases in the LGM wintertime precipitation over the western Mediterranean because of the large-scale effect of the southward shift of the low-level North Atlantic westerly jet associated with changes in the surface air temperature over the northwestern North Atlantic. Some of our LGM experiments show a poleward shift of the North Atlantic wintertime westerly jet, but precipitation increases over the western Mediterranean because of stronger westerlies over the subtropical North Atlantic. These stronger low-level onshore flows, which possibly lessen the severity of winter cold and increase precipitation over western Europe, should be modified.
The Bartlein et al. (2011) data shows the temperatures at 21 ka are similar to or slightly cooler than the 0 ka temperature over Alaska, but our LGM experiments overestimate the annual and wintertime cooling, especially in the GLAC-1D experiment. Ullman et al. (2014) and Liakka and Lofverstrom (2018) found that the LGM Laurentide Ice Sheet elevation increase results in the surface warming over the Arctic, particularly in Alaska and eastern Siberia, because of an increased poleward energy flux by atmospheric stationary waves and radiative effects. Therefore, the GLAC-1D ice sheet, which has the lowest Laurentide Ice Sheet altitude, consistently shows the lowest wintertime temperature in Alaska. However, our overestimation of the Alaskan cooling, in particular the mean annual value, possibly be due to not enough poleward energy flux by atmospheric stationary waves as well as missing some local feedbacks regardless of the height or size of the ice sheets.
Moreover, the reconstruction shows more precipitation in Alaska, while the LGM experiments show decreases in precipitation. The ICE-6G_C experiment shows the precipitation decrease due to shifting the precipitation area to the southern margin of the Laurentide Ice Sheet with winter westerly storm track moving south. This variation has also been reported in previous studies (e.g., Oster et al. 2015). The other two LGM experiments show the decrease in Alaskan precipitation without a southward shift of the westerly storm track. All of our LGM experiments show summer precipitation decrease in Alaska and the increase in western North America due to an equatorward shift of westerlies storm track, and previous studies also show the same pattern (Lofverstrom and Liakka 2016).

Summary
This study comprehensively investigated the impacts of the PMIP4 LGM ice sheets (ICE-6G_C, GLAC-1D, and PMIP3 ice-sheet reconstructions) on the seasonal atmospheric circulation over the Northern Hemisphere using the coupled atmosphere and ocean model HadCM3B-M2.1aD. These ice-sheet reconstructions were created using different approaches, they result in different ice-sheet geometries (and associated boundary conditions). Ice-covered areas of ICE-6G_C and GLAC-1D are similar, but the different distribution of the total ice mass leads to ice sheets with different thicknesses and positions of ice domes, particularly for the North American Ice sheets. The ICE-6G_C and PMIP3 reconstructions have relatively similar ice thickness over the Northern Hemisphere, but the PMIP3 reconstruction has a much more expansive ice sheet over North America. As a result, surface albedo forcing is largest with the PMIP3 ice sheet. The surface albedo feedback is large with the PMIP3 and GLAC-1D reconstructions because of the extended seaice cover over the North Atlantic and North Pacific Oceans. Moreover, the HadCM3B-M2.1aD model simulates colder climate due to a larger sensitivity to changes in the clearsky surface longwave downward radiation than other PMIP3 models.
The differences between the LGM continental ice-sheet geometries have several large and robust impacts on the seasonal atmospheric circulation over the Northern Hemisphere. The Icelandic Low and Azores High in winter and the Pacific High in summer move equatorward. The largescale responses of the seasonal stationary waves to the LGM ice sheets are remarkably similar except for their regional responses. About upper-level (200 hPa) jets, the wintertime North Atlantic het strengthens and extends eastwardly with a larger southwest-northeast tilt in winter, and the summertime North Pacific jet extends across the mid-latitude Pacific Ocean. About low-level (850 hPa) jets, the wintertime North Atlantic jet strengthens and extends across the mid-latitude North Atlantic, and the North Pacific jet strengthens with a larger southwest-northwest tilt and shifts to the equator. Wintertime precipitation increases over western Europe and decreases over storm track regions and the continental ice-sheet areas, while summertime precipitation increases at the west-central USA and decreases over eastern North America and some Northern summer monsoon regions.
Evaluation of our LGM simulations could also provide information about which ice sheet is most realistic. The simulated mean annual temperatures over Europe and Eastern North America at the LGM are adequate; however, the simulated seasonal temperature cycle is inappropriate, especially in Western Europe: underestimating winter cooling and overestimating summer cooling. Compared to the terrestrial-ocean reconstructions, all our LGM experiments tend to overestimate LGM cooling and wet conditions. The LGM experiment with the ICE-6G_C ice sheet simulates the mean climate state closest in amplitude to the mean climate reconstructed from palaeo records, while the LGM experiment with the PMIP3 ice sheet simulates the largest amplitude of this cooling.
The robust patterns of simulated large-scale atmospheric circulations over the Northern Hemisphere can be interpreted as reliable climate responses to the LGM ice sheets. However, the robust responses of seasonal temperature and precipitation in Western Europe, for example, are inconsistent with reconstructions, and our climate model possibly misses some factors such as local feedbacks and hydrological processes over the region. In this way, much paleoclimate research focuses on "testing" the models against palaeodata with an aim to identify problems with the representation of mechanisms within the model (e.g., Flato et al. 2013). However, this is only effective if the magnitude of model-data disagreement exceeds the differences caused by uncertainties in the boundary conditions. Our study suggests that model-data disagreements are robust when considering the global scale (i.e., Had-CM3B-M2.1aD globally overestimates the LGM cooling). However, uncertainties in the ice-sheet area and elevation can significantly alter regional patterns, thus limiting our ability to identify climate model failures at the regional level. Finally, the PMIP4 protocol requires the use of ICE-6G_C, GLAC-1D or the PMIP3 ice sheet for the boundary condition, but considering the above findings, it is much more informative to perform three simulations, one with each ice sheet.