Attribution of the Last Glacial Maximum climate formation

To better understand what determines the Last Glacial Maximum (LGM) cold and dry climate, a suite of numerical experiments with the Nanjing University of Information Science and Technology Earth System Model version 1 are conducted to assess the relative contributions of individual external forcings, including greenhouse gases (GHGs), ice sheets (IS), land-sea configuration (LSC) and the Earth’s orbital parameters, and the contribution of their combination in the LGM mean climate change. Results of the single-forcing sensitivity experiments not only reveal spatial patterns of temperature and precipitation changes that are different from today’s climate, but also shed light on understanding the underlying processes through which each forcing contributes to the formation of the LGM climate. The full forcing experiment simulates a 5.3 K global cooling and an 11.8% reduction of the global mean precipitation, thus yielding a hydrological sensitivity of 2.2% K−1, which is larger than that caused by the present-day GHG forcing. The excessive hydrological sensitivity is primarily attributed to the land-sea configuration change, since its dynamic factor (circulation change) amplifies the precipitation reduction at the tropical convergence zones over both hemispheres. The GHG forcing is the largest contributor to the tropical cooling, whereas ice sheets are responsible for the large hemispheric temperature asymmetry and meridional gradients of the zonal mean temperature change during the LGM period. The LGM precipitation is characterized by decreased precipitation over the Indo-Pacific Ocean and a salient wave train pattern over the Northern Hemisphere (NH). The GHG and LSC forcings are the major contributors to the former since they can change the Indo-Pacific sea surface temperature and the associated Walker circulation, while ice sheets lead to the wave train pattern over the NH by changing the North Atlantic jet stream/storminess and shifting the Intertropical Convergence Zone southward. The climate responses to the LGM forcings are nonlinear. The nonlinearity mainly comes from the overlapping effects induced by the IS and LSC forcings.


Introduction
The Last Glacial Maximum (LGM, approximately 21 ka BP) is the most recent glacial interval during which the global climate is remarkably cold and dry compared to today's climate (Clark et al. 2009). Assessment of the impacts of all external forcings on LGM climatology in coupled climate models has drawn wide attention, especially since the initiation of the Paleoclimate Modeling Intercomparison Project (PMIP, Joussaume and Taylor 1995). This project 1 3 could be characterized by the atmospheric general circulation models (AGCMs) that are forced by the prescribed SST and sea ice concentration in PMIP phase 1 (PMIP1), coupled ocean-atmosphere/ocean-atmosphere-vegetation GCMs in PMIP2 (Braconnot et al. 2007a, b), and the Earth System Models (ESMs) in PMIP3 (Harrison et al. 2015).
It is worth noting that the afore-mentioned studies have been carried out in a multi-forcing framework without isolating impacts of individual climate forcings. However, understanding the contribution of each single forcing to the LGM climate formation is more useful for understanding climate dynamics. Technically, single forcing experiment is a tractable way in climate modeling (Kageyama et al. 2017).
Separate investigation of the impact of each forcing on the global scale can be traced back to Manabe and Broccoli (1985) and Broccoli and Manabe (1987). They used a simplified AGCM and coupled it with a slab ocean model to examine the role of ice sheets and GHGs in LGM climate formation. They suggested that the GHG and ice sheet forcings play equally important roles in the global mean surface temperature changes. However, lacks of ocean dynamics and cloud feedback could lead to underestimation of the amplitudes of the responses (Webb et al. 1997;Justino et al. 2005). In coupled GCMs, it was found that the lower CO 2 forcing during the LGM accounted for about half of the temperature decrease, while the effects of ice sheets and land-sea configuration were only inferred from the difference between the LGM full forcing experiment and the lower CO 2 simulation (Kim 2004;Laîné et al. 2009a;Brady et al. 2013).
On regional scales, Kim et al. (2017), using an AGCM to isolate the impacts of individual boundary conditions, found that the sea ice forcing in the Southern Ocean is the most important factor that controls the strength and location of the Southern Hemisphere (SH) westerlies. Over the NH, different ice sheet topographies were specified in a high-resolution AGCM and the importance of the Laurentide ice sheet was identified in altering the North Atlantic jet stream and the precipitation pattern over North America and Europe (Pausata et al. 2011;Hofer et al. 2012). Di Nezio et al. (2016) investigated the rainfall pattern over the Indo-Pacific warm pool with the Community Earth System Model version 1 (CESM1) by altering the Maritime Continent (MC) land-sea configuration. Their results suggested that continental shelf exposure plays a critical role in explaining LGM precipitation anomalies over the Indian Ocean.
Most of the previous works examining the influence of each single LGM forcing are focused on regional scales and use AGCMs alone or simplified coupled GCMs. They do not identify the forcings' impacts in a coupled GCM framework and on a global scale. Single forcing experiment with coupled climate models is advantageous for assessing the relative contributions of each forcing to global climate change, which includes changes in the mean states of precipitation and temperature, global monsoon, ITCZ and ENSO. For example, Yan et al. (2016Yan et al. ( , 2018 found that the Australian monsoon, in contrast to all other regional monsoons, is strengthened and expanded during the LGM. This unique feature cannot be qualitatively explained by a full forcing experiment or the GHG forcing alone. Identification of key factors responsible for this unique feature relies on understanding each single forcing's contributions on the global scale in coupled GCMs. Single forcing experiment with global model can also facilitate identification of the major process responsible for hydrological sensitivity. It is the slope of globally averaged precipitation change with respect to surface temperature change that reflects the global hydrological cycle change in response to temperature change (eg. Held and Soden 2006;Bala et al. 2008). Recent studies reported that hydrological sensitivity during the LGM is 15-23% higher than the GHG-induced hydrological change in PMIP3 models (Li et al. 2013;Harrison et al. 2015;Yan et al. 2016). However, it remains unclear what causes the larger hydrological sensitivity during the LGM. According to the Clausius-Clapeyron equation, temperature increase would raise the saturation vapor pressure at a rate of about 7% in the current climate, but the climate model-simulated hydrological sensitivity is only about 1-2% K −1 , implying that the change of global circulation should account for the rest of the change (about 5% K −1 ) (Held and Soden 2006). From an energetic consideration, this sensitivity may depend on the mean state (Allan 2009;Allen and Ingram 2002), and may also be affected 1 3 by different external forcings (Kleidon and Renner 2013). Therefore, it is of interest to figure out which external forcing causes the larger hydrological sensitivity during the LGM.
In this study, a suite of single-and combined-forcing sensitivity experiments are designed and conducted in order to be compared with the preindustrial (PI) control experiment. The main purpose of this study is to access the impacts of each individual forcing on the change of LGM temperature and precipitation patterns and the hydrological sensitivity. This paper is organized as follows. The model, experimental design, and model validation are described in Sect. 2. The attribution of LGM temperature formation is presented in Sect. 3. The characteristics of precipitation responses to different forcings and the hydrological sensitivity are analyzed in Sects. 4 and 5, respectively. Major findings are summarized in Sect. 6.

Model and experimental design
The model used is the recently developed Nanjing University of Information Science and Technology Earth System Model version 1 (NESM v1, Cao et al. 2015). The model is designed to study interactions between different components of the earth climate system and its responses to natural and anthropogenic forcings. NESM v1 consists of the atmospheric, land, ocean and sea-ice components. The atmospheric component is the European Centre Hamburg Model (ECHAM v5.3) (Roeckner et al. 2003). It is an AGCM with a spectral dynamical core formulated in spherical harmonics. The atmospheric resolution used here is T42L31, which has a horizontal resolution of ~ 2.8° latitude by 2.8° longitude and 31 vertical levels extending from the surface to 10 hPa. The oceanic component is the Nucleus for European Modeling of the Ocean (NEMO), which represents threedimensional circulation of the global ocean. NESM employs OPA (Ocean PArallelise) of the NEMO framework for the ocean dynamics and thermodynamics (Madec and NEMO Team 2012). NEMO runs with an ORCA2 grid which has 31 vertical layers with 9 layers in the upper most 100 m and an inhomogeneous horizontal resolution with higher resolution in low latitude regions (0.5°-2°). The sea-ice component of NESM is the Los Alamos sea-ice model (CICE v4.1), which includes four vertical ice layers and one snow layer and adopts a multi-layer thermodynamic scheme (Hunke and Lipsomb 2010). Performances of NESM in simulating modern climate were evaluated under the modern climatological forcing. It reproduces reasonably well observed large-scale present-day climatic fields such as SST and precipitation (Cao et al. 2015. NESM also captures observed larger-scale climate variability and climate modes such as El Niño-Southern Oscillation (ENSO) and Madden-Julian oscillation (MJO) reasonably well.
The LGM external forcings are significantly different from the current ones. Four types of forcings are considered in this study following the PMIP3 protocol, including the Earth's orbital parameters, greenhouse gases, ice sheets, and land-sea configuration. Besides, the runoff, land surface vegetation and aerosol are specified as the PI condition. The additional continental shelves due to the lower sea level are set to have the same vegetation types as their neighboring grid cells in the zonal direction. The ocean salinity is recommended to be added by 1 psu, while it is not considered in the simulations. Change of the Earth's orbital parameters affects the seasonal cycle of solar insolation, but the change in global mean insolation is negligibly small (Berger 1978). Decreased GHG concentration would reduce net radiative forcing, and the presence of continental ice sheets is associated with the lowering of sea level, altering of land-sea configuration and land surface elevation (Fig. 1). The LGM ice sheet forcing is a blended product obtained by averaging three different ice sheet reconstructions: ICE-6G v2.0 (Argus and Peltier 2010), Meltwater routing and Ocean-Cryosphere-Atmosphere response (MOCA; Tarasov and Peltier 2004), and Australian National University (ANU; Lambeck et al. 2010). Ice sheets are treated as mountain ranges with glacial surface properties (e.g. albedo). The mean sea level drop is about 120 m in the reconstruction data of ice sheets (Brady et al. 2013;Kageyama et al. 2017). The land area is increased by about 5% at the latitudes, except over 70°-80°N where the land area is increased by more than 50% (Fig. 1b). The changes of land-sea configuration mainly include the exposure of the Sunda and Sahul Shelf, and the closure of Bering Strait.
To quantitatively estimate the climate effects of external forcings, a preindustrial (PI) control experiment and five sensitivity experiments as well as the LGM full forcing experiment are designed. The boundary conditions and forcings of the PI and LGM experiments follow the PMIP3 experimental protocol. Four single-forcing sensitivity experiments, which are designed based on the PI control experiment, are named as the Earth orbit (EO), greenhouse gas (GHG), ice sheet (IS), and land-sea configuration (LSC) experiments. In the EO experiment, the Earth's orbital parameters, including eccentricity, obliquity, and perihelion, are specified as those in the LGM conditions. In the GHG sensitivity experiment, the concentrations of CO 2 , CH 4 and NO 2 are lowered to 185 ppm, 350 ppb and 200 ppb, respectively. In the IS experiment, the LGM continental ice sheet cover is used. The changes of ocean bathymetry and land surface elevation are only considered over the specified ice sheet regions. The 1 3 land sea configuration is specified to the LGM state in the LSC experiment, while the continental ice sheets are the same as those in the PI experiment. This is targeted at investigating the influence of exposed land shelves. Although the LSC and IS experiments are independently considered, since the presence of ice sheets, lowering of sea level and change of land-sea configuration are intrinsically related, it is meaningful to consider their combined effects. Therefore, an additional sensitivity experiment is carried out in which the combined effects of the ice sheet and land-sea configuration forcings are considered together. It is named as the ice sheet-land-sea configuration (IS_LSC) experiment. Detailed experimental design and spin-up length of all simulations are summarized in Table 1.
The initial conditions for the PI control run are extracted from an equilibrium state of NESM v1 which is forced by the external forcings fixed at the level of the year 1990 (Cao et al. 2015). Before the PI experiment, a 500-year coupled model spin-up is conducted with the constant PI forcing. The initial conditions of the LGM and all sensitivity experiments are derived from the PI experiment. Figure 2 shows the time evolution of the global mean surface air temperature (TAS) and subsurface ocean (0-600 m) temperature in all simulations. To obtain a quasi-equilibrium state for each simulation, we use the linear trend of the globally averaged TAS as the matric, although the global mean SST is a popular metric for model stability in previous paleoclimate studies LGM ice sheets, and the shaded area is ice sheet topography. b The zonal mean fraction of expanded land surface area 1 3 (Braconnot et al. 2007a;Zhang et al. 2013). Zhang et al. (2013) indicated that the surface characteristics can be similar and quasi-stationary, even when a significant trend in deep ocean properties still exists. Indeed, the subsurface ocean (0-600 m) temperature is still drifting in LSC_IS, IS and LSC experiment, although the spin-up length exceeding 700-year. In PI, LGM, EO and GHG experiment, the subsurface ocean is relative stable. In all experiment, the TAS is in a quasi-equilibrium state. Here, the quasi-equilibrium state requires the global mean TAS trend to be less than 0.02 K century −1 . This criterion is comparable with the temperature linear trend in the PI experiments of most of the Coupled Model Intercomparison Project Phase 5 (CMIP5) models (Gupta et al. 2013). In this study, the last 100-year data from each simulation are analyzed.

Model validation
The TAS and precipitation are the most important fields in the coupled system. Figure 3 shows the climatological mean TAS and precipitation simulated in the PI and LGM experiments as well as their differences. The coupled model realistically simulates the annual mean temperature of the PI experiment ( Fig. 3a), including the Indo-Pacific warm pool and the eastern Pacific cold tongue, as well as the temperature in the plateau regions with high elevations. The simulated LGM climate is colder and drier than the PI climate (Fig. 3e, f). The omission of vegetation change and aerosol change in the experimental design could induce underestimation of the simulated surface cooling, since the negative radiative feedback could further cool the Earth's surface (Crucifix and Hewitt 2005;Schneider von Deimling et al. 2006a, b). A salient feature is stronger cooling over land than over ocean and amplified cooling over high latitudes, especially over the NH land. Over the high latitudes, the temperature decreases by more than 12 K due to the presence of continental ice sheets and ice-albedo feedback induced "polar amplification". The greatest temperature decreases are located over the regions occupied by the Fennoscandian and Laurentide ice sheets with the maximum decrease reaching up to 18 K over the North American continent. It is the result of elevated topography and increased surface albedo. Interestingly, there is a warm signal over the Davis Strait and Labrador Sea, which is also seen in LGM proxy data (MARGO Project Members 2009; Annan and Hargreaves 2013) and simulated by CMIP5 models (e.g., IPSL-CM5A, Kageyama et al. 2013). The sea ice transport cut-off over the Northern Canada and the enhanced North Atlantic subtropical high both contribute to the warmer SAT. The model also captures the major features of the spatial distribution of climatological precipitation in the PI experiment, including precipitation over the oceanic convergence zones, mid-latitude storm tracks and the less precipitation regions in the subtropics (Fig. 3b). In the LGM experiment, the change of precipitation shows a clear "wet-gets-dry" pattern (Held and Soden 2006;Wang et al. 2012) and an overall decrease on the global scale (Fig. 3f). Precipitation is significantly reduced over the ITCZ and South Pacific Convergence Zone (SPCZ) as well as over the mid-latitude storm track regions where the moisture convergence and upward motion are strong. Over the subtropics, the moisture divergence is reduced and thus precipitation is increased, although the underlying SST is colder than that in the PI control run, suggesting the dynamic control over precipitation change. In addition, in the LGM experiment the simulated precipitation is reduced over the western European and northern American continents, where low-level anticyclonic circulation anomalies dominate (Fig. 3f) due to the presence of ice sheets. In contrast to the ice sheet regions, the simulated precipitation is increased to the south of the Fennoscandian and Laurentide ice sheets. Contrary to the high-latitude dry climate in the LGM experiment, moisture convergence and a wet climate can be found over the Labrador Sea due to an anomalous cyclonic circulation. Overall, the simulated LGM precipitation pattern is consistent with that in the PMIP3 models and the LGM proxy data (Brady et al. 2013;DiNezio et al. 2011;Oster et al. 2015).

Surface temperature responses to LGM forcings
The surface temperature responses are examined in all sensitivity experiments. Figure 4 shows the temperature changes relative to the PI experiment in the LGM, EO, GHG, IS_ LSC, IS and LSC experiments, respectively. Temperature changes induced by the Earth's orbital parameter changes are negligibly small compared to the responses to other forcings (Fig. 4b). This is because the global mean insolation increases only by about 0.03 W m −2 , although there is an apparent insolation change in the seasonal cycle (Fig. 5a). The incoming solar radiation is decreased by more than 8 W m −2 over the high latitudes during the local summer, leading to the larger temperature responses over the high latitudes of both hemispheres (Fig. 5b). The high latitude temperature is characterized by the local summer cooling and local winter warming. In addition, the temperature response is delayed to insolation change. The decrease of GHG concentrations significantly affects the global surface temperature (Fig. 4c), which is consistent with a net radiative forcing decrease of 2.8 W m −2 induced by the lower GHG forcing. At the same latitude, the response over the land surface is stronger than over the ocean. The TAS reduction exceeds 5 K in the entire Arctic Ocean, while large reduction of temperature is only seen in part of the Ross Sea and Weddell Sea over the Antarctic, indicating a stronger polar amplification in the NH than in the SH. In contrast, temperature change induced by the IS_LSC forcing is much larger in high latitudes than in the tropics; it also has a more complex regional structure in high latitudes (Fig. 4d), especially over NH high latitudes. Salient temperature decreases are found over North America, Northern Europe, Antarctica and the Kuroshio extension region. The surface temperature is increased over the Davis Strait, Labrador Sea and Bering Strait. This pattern of temperature change is also presented in the LGM full forcing experiment (Fig. 4a), suggesting a critical role of the IS_LSC forcing in formation of the NH high latitude temperature pattern. The IS_LSC forcing can be further decomposed to the ice sheet (IS) and land-sea configuration (LSC) forcings. In the IS experiment, the strongest temperature responses mainly appear over mid-and high latitudes (Fig. 4e), which coincide with the Fennoscandian and Laurentide ice sheets over the European and American continents, the continental ice topography over the Antarctic continent and the Patagonian ice sheet over the Southern American continent. Meanwhile, temperature decreases more in the NH than in the SH. The temperature pattern also suggests that the warm signal over the Labrador Sea and Bering Strait in the IS_LSC experiment could be partially traced back to the IS forcing. The tropical cooling in the LSC experiment resembles that in the IS_LSC experiment, which means that the LSC forcing contributes to the tropical cooling. As expected, temperature decrease due to the LSC forcing (Fig. 4f) exceeds 1 K over the tropics, and the maximum cooling occurs where the land-sea configuration has dramatic change. The TAS decreases by more than 2 K around Antarctica, the MC, western Europe, the Canadian Arctic Archipelago and the coast of East Asia. In contrast, the TAS increases by about 5 K over the Arctic (Fig. 1b). The change of surface status as well as albedo would increase the absorbed solar radiation and warm the Arctic. In summary, influence of the GHG forcing is relatively uniform over the globe compared to that of the other forcings, and the special regional patterns of temperature change in the LGM are largely attributed to the IS forcing.
The zonal mean temperature changes in each experiment are shown in Fig. 6. Between 60°S and 30°N, the temperature decreases by about 3.5 K in the LGM experiment (Fig. 6a, black line), which can be attributed to the GHG (Fig. 6a, purple line) and IS_LSC forcings (Fig. 6a, red line) as the EO forcing (Fig. 6a, yellow line) has nearly imperceptible impacts over all latitudes. The temperature reductions gradually increase poleward starting from 30°N and 60°S. In 1 3 the SH high latitudes (60°S-90°S), the latitudinal distribution of LGM temperature change is dominated by the IS_ LSC forcing. Over the NH high-latitude, the IS forcing is the dominant contributor to the zonal distribution of temperature change (Fig. 6b), although it is partially compensated by the impact of the LSC forcing. This suggests that the IS forcing is effective in altering the meridional temperature gradients. Table 2 summarizes the globally and hemispherically averaged temperature changes relative to the PI control run in the LGM and five sensitivity experiments. The global mean temperature changes by 0.1 K, − 2.8 K and − 2.5 K in the EO, GHG, and IS_LSC experiments, respectively. The sum of global mean values is close to that in the LGM full forcing experiment. Examining the IS and LSC forcings separately, however, we found that the temperature decreases by 1.2 K and 2.5 K respectively, indicating an overestimation of the reduction due to the overlapping (nonlinearity) effects that occur when the IS and LSC forcings are simultaneously imposed. On the hemispheric scale, the NH is 1.7 K cooler than the SH in the LGM experiment. This is due to the difference in polar amplification, as it extends poleward from 30°N in the NH but from 60°S in the SH, respectively. The hemispheric temperature asymmetry is more significant in the IS_LSC experiment than in the GHG experiment. The main contributor is the IS forcing which produces a hemispheric temperature asymmetry of 1.4 K, although it is partially offset by the LSC forcing ( Fig. 6; Table 2). Figure 7 shows the changes of 850 hPa winds and precipitation in the LGM full forcing experiment and each sensitivity experiments with reference to the PI control run. In the LGM experiment (Fig. 7a), precipitation reduction is most prominent over the MC and the adjacent Indo-Pacific Ocean warm pool, as well as over the NH ice sheet regions. This is consistent with the results from PMIP2 and PMIP3 multi-model ensemble mean (  ice sheet and land-sea configuration (IS_LSC) forcing is the primary contributor, except in the equatorial western Pacific where the GHG forcing also plays an important role. Consistent with its role in temperature change, the EO forcing is a negligible contributor to circulation and precipitation changes (Fig. 7b). Different single-forcing experiments exhibit distinguished precipitation change patterns. In the GHG experiment, the  . 7 Spatial maps of mean annual precipitation (mm day −1 ) and 850 hPa winds (m s −1 ) differences relative to the PI control run for the LGM (a), EO (b), GHG (c), IS_LSC (d), IS (e) and LSC (f) experiments, respectively. Only value that are statistically significant at 5% confidence level based on the two-sided Student's t-test are shown. The grids with elevations higher than 1500 m are ignored 1 3 simulated precipitation changes show a robust "wet-getsdry" pattern in the cold LGM climate (Fig. 7c). The most pronounced precipitation reduction occurs over the equatorial western Pacific Ocean, although precipitation is also decreased over mid-latitude storm track regions. Positive precipitation anomalies are shown in subtropical divergent regions. This precipitation change pattern resembles the one projected by the CMIP5 model but they have opposite polarities (Lee and Wang 2014).

Attribution of precipitation change to individual LGM forcings
In the IS_LSC experiment, on the other hand, the precipitation anomaly pattern resembles the result from the LGM experiment except over the equatorial western Pacific Ocean (Fig. 7d). The IS_LSC forcing produces a NH mid-latitude wave-train pattern and a tropical east-west sandwich precipitation anomaly pattern (decreased precipitation over the MC and its vicinity and increased precipitation over the western Indian Ocean and western equatorial Pacific). These two precipitation patterns are shown in the IS and LSC experiments respectively, indicating that different physical processes are involved in forming different regional precipitation anomalies under the IS forcing and LSC forcing. And the two forcings are largely responsible for the NH mid-latitude wave train pattern and the tropical east-west sandwich precipitation pattern during the LGM period.
The IS forcing generates anticyclonic circulation anomalies over the two ice sheet regions and the subtropical Atlantic Ocean, and produces cyclonic circulation anomalies over the Northern Pacific and the east coast of North America (Fig. 7e). The wave train circulation anomaly pattern resembles the response to the orographic forcing of an idealized entropic mountain (Hoskins and Karoly 1981). Furthermore, complex GCM studies also confirm that the ice sheet topography accounts for a larger portion of circulation changes over the midlatitudes, indicating the important role of the orographic effect of ice sheets (Chiang et al. 2003). The precipitation anomalies correspond to the circulation anomalies well with the anticyclonic anomalies and cyclonic anomalies. Moreover, the IS forcing reduces the ITCZ precipitation, specially over the western hemisphere (Fig. 7e). In the LSC experiment, precipitation response is most pronounced over the equatorial region (Fig. 7f) due to the local and/or remote response to exposure of the Sunda Shelf and Sahul Shelves. The east-west sandwich pattern of precipitation anomalies, which is clearly shown in this experiment, is associated with the anomalous equatorial Indian Ocean easterlies and Pacific westerlies. Meanwhile, the cyclonic vorticity associated with equatorial westerlies over the Indian Ocean enhances Indian monsoon rainfall.
In sum, the ice sheet forcing is a dominant factor in shaping the NH high latitude precipitation change, while the tropical precipitation response results from the combined impacts of the land-sea configuration and greenhouse gas concentration forcings, given the negligibly small influence of the EO forcing.

Forcing mechanisms responsible for LGM precipitation changes
To understand the mechanisms behind precipitation responses to the GHG and LSC forcings, we investigated the related atmospheric lower boundary condition changes because tropical precipitation is sensitive to the underlying SST gradient. The tropical SST changes in the GHG and LSC experiments are compared in Fig. 8. In the GHG experiment (Fig. 8a), the simulated SST shows a La Niñalike cooling with large temperature gradient over the central Pacific. The La Niña-like cooling enhances easterlies over the Pacific and westerlies over the northern equatorial Indian Ocean. Therefore, the ascending branch of the Walker circulation is enhanced to the west of 130°E and strong Fig. 8 The changes of tropic SST (K) and 850 hPa winds (m s −1 ) in GHG and LSC experiments relative to PI control run. Only value that are statistically significant at 5% confidence level based on the two-sided Student's t-test are shown 1 3 subsidence anomaly occurs over the western Pacific (Fig. 9b) where precipitation is severely suppressed (Fig. 7c). Over the MC, the effects of strengthened ascending motion are partially compensated by the decreased boundary layer moisture content that is induced by a cold SST condition, yielding a modest precipitation change. However, over the equatorial Western Pacific (EWP), the anomalous subsidence and cold SST work together, resulting in a dramatic precipitation reduction (Fig. 7c), which is consistent with the analysis of Chadwick et al. (2013). This mechanism is responsible for the precipitation reduction over the EWP in the GHG experiment. When the LGM land-sea configuration is imposed, the exposure of land shelves substantially reduces surface evaporation and moisture supply to precipitation, thus dramatically reducing precipitation over the expanded land (Fig. 7f). The loss of precipitation heating over the expanded land produces equatorial westerly anomalies over the Indian Ocean due to a descending Rossby wave response and easterly anomalies over the equatorial western Pacific which can be viewed as a Kelvin wave response (Fig. 8b). The precipitation and circulation anomalies caused by the LSC forcing tend to have an opposite sign from those caused by the GHG forcing over the Indian Ocean and western Pacific (Figs. 7c,  f, 8a, b). The easterly anomaly over the eastern Indian Ocean may favor coastal upwelling along the west cost of Sumatra and eastern equatorial Indian Ocean, leading to more SST cooling over the eastern Indian Ocean. However, over the western Indian Ocean SST is less affected (Fig. 8b). The resultant westward SST gradients would further enhance the easterly wind anomaly, so that the SST gradients would be amplified and maintained by the Bjerknes-type of positive feedback (Bjerknes 1969). Dynamically, the anomalous Indian Ocean SST gradient accelerates the low-level divergence flow, and the boundary layer divergence weakens the ascending branch of the Walker circulation over the MC (Fig. 9c). Consequently, precipitation is reduced over the equatorial Indian Ocean and the adjacent land. Over the EWP, the anomalous westerly wind enhances the low-level convergence and ascending motion (Fig. 9c), therefore increasing precipitation. However, the precipitation response is much weaker compared to the GHG-induced precipitation reduction; thus, the EWP precipitation change is ruled by the GHG forcing. In summary, it is the combined effect of the GHG and LSC forcings that causes anomalous Walker circulation and precipitation reduction over the whole Indo-Pacific region during the LGM period.
Over the NH mid-latitudes, the atmospheric circulation and hydrological cycle are substantially altered by the imposed glacial boundary conditions, which creates a wave-like circulation pattern with strong anticyclonic circulations sitting on the ice sheets, and shifts winter Aleutian Low southward (figure not shown). Significant decrease of surface temperature and the anomalous anticyclonic circulations induced by low-level divergence reduce surface evaporation and low-level moisture convergence, therefore reducing local precipitation. On the other hand, the ice sheet forcing-related changes in the semi-permanent system are likely to be responsible for the eddy-driven jet stream change (Laîné et al. 2009b;Hofer et al. 2012;Oster et al. 2015;Wang et al. 2018). Figure 10 shows the NH winter (DJF) 850 hPa wind difference between the IS experiment and the PI experiment. Significant wind changes are located in the North Pacific and the North Atlantic. The westerlies are enhanced/weakened over the south/north of the westerly jet over the North Pacific sector, indicating the southward shift of the North Pacific Fig. 9 Changes in Walker circulation (shaded, hPa day −1 ) in response to LGM, GHG and LSC forcing. The contours outline the mean Walker circulation in PI control run. Only value that are statistically significant at 5% confidence level based on the two-sided Student's t-test are shown 1 3 jet stream. Correspondingly, accelerated zonal wind extends from the south of the Laurentide Ice Sheet to the Mediterranean region. It enhances the North Atlantic jet stream and shifts it southward. The equatorward shift of the westerly jet in the IS experiment is similar to that in the full LGM experiment (figure not shown) and the PMIP3 MME results (Oster et al. 2015;Wang et al. 2018), suggesting the central role of continental ice sheets in modulating the NH westerlies. The presence of ice sheets alters the meridional temperature gradient which is closely related to the eddy baroclinic activity. Precipitation associated with the eddy baroclinic activity is enhanced over the south of ice sheets . With the eastward extension of the westerly jet, the storm activity is increased (decreased) over the southern (northern) side of the Fennoscandian ice sheet (Hofer et al. 2012), leading to more vigorous ascending motion in southern Europe and enhanced subsidence in northern Europe. This corresponds to the increased (decreased) precipitation over the southern (northern) Europe. Figure 11 shows the zonal mean precipitation changes simulated in the LGM full forcing experiment and five sensitivity experiments. In all but the EO experiment, decreased precipitation occurs over the equatorial and storm track regions, where the climatological zonal mean precipitation reaches maxima. Moreover, the reduction of precipitation in the LGM experiment is more evident over the equatorial NH than over the equatorial SH. The contributions of each forcing to zonal mean precipitation are different. The GHGinduced precipitation change is more concentrated at the equator (purple line), while the precipitation change caused by the IS_LSC forcing is located over the convergence zones of both hemispheres (red line). More specifically, the zonal mean precipitation response in the IS_LSC experiment arises from (1) reduction of the NH ITCZ precipitation due to the IS forcing (green line) and (2) reduction of convergence zone precipitation in both hemispheres due to the LSC forcing (blue line).
To understand the physical processes responsible for the precipitation change, we investigated the column-integrated moisture tendency equation (Hsu et al. 2012).
where ⟨q⟩ is the total column-integrated water vapor, ∂/∂t is the time tendency, 〈 〉 indicates the vertical integration from the surface to 100 h Pa, ∇ is the gradient operator, q is the specific humidity, V is the wind vectors, and E and P are the surface evaporation and precipitation. Since all simulations are in a quasi-equilibrium state, the term ∂w/∂t Fig. 10 a The Northern Hemisphere winter (DJF) 850-hPa winds (m s −1 ) and zonal wind speed (shaded) in PI experiment, and b the difference between IS experiment and PI experiment, only value that are statistically significant at 5% confidence level based on the two-sided Student's t-test are shown. The grids with elevations higher than 1500 m are ignored. The purple lines outline the LGM ice sheets Fig. 11 Zonal mean precipitation anomalies relative to PI control run for a LGM, EO, GHG and IS_LSC experiment, and b IS_LSC, IS and LSC experiment. The solid lines indicate that the values statistically significant at 5% confidence level based on the two-sided Student's t-test, while the insignificant values are indicated by dashed lines 1 3 would vanish. The change of monsoon precipitation could be further derived as: where ∆ means the difference between two experiments. Again, 〈 〉 indicates the vertical integration from the surface to 100 h Pa. The first term in the right-hand side is moisture advection, and the second and third terms indicate the contribution of moisture convergence and surface evaporation, respectively. We ignored the moisture advection term due to its small contribution and the evaporation term because it is meridionally uniform over tropics. The changes of moisture convergence could be further decomposed to the circulation change-related dynamic process, temperature change-related thermodynamic process, and their interactions. They are derived as follows: where ∆ means the difference between the forced experiment and the PI control experiment, the variable with subscript 'pi' denotes the PI experiment, and D indicates the divergence. The first term at the right-hand side of Eq. (3) is associated with the circulation change, which could be regarded as the dynamic contributor. The second term reflects the change of moisture, which could be regarded as the thermodynamic effect. Figure 12 shows the zonal mean thermodynamic and dynamic effects as well as precipitation changes in each experiment, and the red dots indicate the tropical mean evaporation. Since the surface evaporation is zonal uniform, the meridional distribution of precipitation ( P ) would be consistent with the meridional distribution net precipitation ( P net ). The sums of thermodynamic and dynamic effects (solid black lines) coincide well with the meridional distribution of precipitation anomalies (dashed line), especially over the tropics. All experiments show similar thermodynamic effects (blue lines), whereas the dynamic effects (red lines) vary among different experiments. Enhanced convergence produces increased equatorial precipitation under the EO forcing. In the GHG experiment (Fig. 12c), the influence of the vertical motion confirms the "wet-gets-dry" mechanism, which reproduces the enhanced vertical motion in the NH ITCZ (Fig. 13a). It compensates the reduction due to decreased boundary layer moisture content, yielding a modest precipitation change near 12°N. However, the ascending motion is suppressed near the equator (5°S-5°N), which reduces the equatorial precipitation in the GHG experiment (Fig. 13a). The dynamic effect of the IS forcing weakens the NH ITCZ but strengthens the SH convergence zone (Fig. 13b), which is induced by the hemispheric temperature asymmetry ( Table 2). The NH cools more than the SH does under the IS forcing, causing southward cross-equatorial flows that enhance convergence and precipitation over the SH convergence zone (Fig. 12d). This leads to the southward shift of Hadley circulation (Fig. 13b). Under the LSC forcing, the surface condition change induced air-sea interaction over the equatorial Indian Ocean suppresses the local ascending motion (Fig. 12f), resulting in a large precipitation decrease over the SH convergence zone. A weakened NH ITCZ is also shown, whose cause needs further investigation. In summary, the sum of dynamic and thermodynamic effects could well explain the meridional distribution of precipitation changes in the sensitivity experiments. The different locations of maximum precipitation change among sensitive experiments are due to the variant dynamic processes induced by different forcings. Table 3 summarizes the global and hemispheric mean precipitation changes. In the LGM experiment, the global mean precipitation decreases by 0.34 mm day −1 , which leads to a global mean hydrological sensitivity of 2.2% K −1 (Table 4). This result is consistent with the PMIP3 model results (Li et al. 2013;Jiang et al. 2015;Yan et al. 2016). Li et al. (2013) reported the range of values obtained from PMIP3 model is between 1.80 and 2.89% K −1 . The multi-model ensemble mean value is 2.26% K −1 (Jiang et al. 2015;Yan et al. 2016). However, the future projected hydrological sensitivity is smaller in CMIP 5 models.

Hydrological sensitivity
The decrease in global mean precipitation during the LGM is close to the sum of reductions induced by the GHG and IS_LSC forcings (Table 3). About 43% of the LGM precipitation change can be attributed to the change of GHG concentration, while the IS_LSC accounts for about 57%. Precipitation decrease under the IS_LSC forcing (− 0.21 mm day −1 ) is larger than that under the GHG forcing (− 0.16 mm day −1 ). However, temperature changes induced by the GHG and IS_LSC forcings are almost the same. Therefore, the IS_LSC forcing induces larger hydrological sensitivity than the GHG forcing does. In the LGM experiment, the reduction of precipitation in the NH is 50% higher than that in the SH (Table 3), which arises from the fact that the IS forcing-induced precipitation change occurs mainly over the NH (Table 3).
Precipitation increases by 3.9% when solar radiation warms the Earth's surface by 1 K in the EO experiment (Table 4). This is close to the estimate made in previous studies (Tilmes et al. 2013;Kleidon and Renner 2013). Given the small contribution of the EO forcing, the high hydrological sensitivity in the LGM experiment comes primarily from the IS_LSC forcing. Comparison of the results from single forcing experiments with those in the IS_LSC forcing experiment reveals that the LSC forcing is the major 1 3 contributor to the high hydrological sensitivity during the LGM. The simulated hydrological sensitivity is about 4% K −1 in the LSC experiment.
The hydrological sensitivity varies with latitude and large sensitivity is seen over the tropics. This is because the change of precipitation is large in the latitudes with heavy precipitation (tropics) (Fig. 11). Yet the temperature change is small over the tropics (Fig. 6). This assertion is supported by the zonal mean hydrological sensitivities as shown in Fig. 13. As expected, the tropical region has large sensitivities, especially in the IS_LSC experiment (Fig. 13a). It is also illustrated by the result of the IS and LSC experiments (Fig. 13b). The IS forcing amplifies the hydrological sensitivity over the NH ITCZ, but reduces the sensitivity in the SH convergence zone, resulting in a modest global mean sensitivity (Table 4). However, the LSC forcing produces large hydrological sensitivity in the convergence zones of both hemispheres, yielding a high global mean value (Table 4).
Compared to the GHG forcing, the IS and LSC forcings generate higher hydrological sensitivity over the tropical convergence zones in both hemispheres. This could be Fig. 12 Zonal mean precipitation (dashed lines, mm day −1 ), moisture convergence (black lines, mm day −1 ) that induced by moisture anomaly (blue lines, mm day −1 ) and induced by circulation anomaly (red lines, mm day −1 ) relative to PI control run for LGM (a), EO (b), GHG (c), IS_LSC (d), IS (e) and LSC (f) experiment, respectively. The red dots on the right y-axis indicate the global mean surface evaporation 1 3 attributed to the differences in the meridional distribution of dynamic processes, which are closely linked with locations of the effective external forcing, as discussed in the previous section (Fig. 14).

Summary
Assessment of the influences of the GHG concentration (GHG), ice sheet (IS), land-sea configuration (LSC) and the Earth's orbital (EO) forcings is imperative for understanding the LGM climate formation. In this study, the PI control run, the LGM full forcing experiment and five sensitivity experiments are conducted with the NESM v1 model. Results of the PI control run and the LGM experiment demonstrate that the NESM v1 model can reproduce reasonable climate mean states and responses to large external forcings, especially in the temperature, precipitation and circulation fields. The LGM full forcing experiment estimates that the global mean surface temperature cooled by 5.3 K, and the global mean precipitation rate decreased by 0.34 mm day −1 , resulting in a global precipitation decrease at a rate of 2.2% K −1 . This is consistent with previous proxy data and results of numerical simulation studies. For single-forcing experiments both the sum of temperature and that of precipitation are larger than their counterparts in the LGM full forcing experiment. Further analysis shows that the nonlinearities come from the overlapping effects induced by the IS and LSC forcings (Tables 2, 3).
Results of the single-forcing sensitivity experiments reveal the spatial patterns of temperature and precipitation changes induced by each external forcing. In terms of the temperature change, the lowered GHG concentration induces a pattern that is similar to the global warming pattern but has opposite polarities. The presence of ice sheets causes the local and downstream cooling, and it is responsible for the large northern and southern hemispheric temperature asymmetry and the meridional gradients of the zonal mean temperature changes during the LGM. The alteration of land-sea configuration cools the expanded land areas by changing the surface albedo, but could warm the NH high latitudes due to the lower boundary that changes from sea ice to land surface.

Fig. 13
Changes in Walker circulation (shaded, hPa day −1 ) in response to a GHG and b LSC forcing. The vectors are the composite of zonal wind (m s −1 ) and vertical motion (ω × 100, unit: hPa day −1 ). Only value that are statistically significant at 5% confidence level based on the two-sided Student's t-test are shown Table 3 The global and hemispheric averaged precipitation (mm day −1 ) differences relative to the PI control run for the LGM and five sensitivity experiments The values statistically significant at 5% confidence level based on the two-sided Student's t-test are bolded

3
In terms of precipitation change, the LGM precipitation response is characterized by a dry Indo-Pacific Ocean and a salient wave train pattern over the NH. The former is mainly caused by the GHG and LSC forcings, while the NH wave train pattern is forced by the land ice sheets. The GHG forcing significantly reduces the EWP precipitation, and the exposure of land shelves weakens precipitation over the MC region. This suggests that the Indo-Pacific warm pool precipitation anomaly in the LGM is due to the combined effects of the LSC and GHG forcings. The presence of Fennoscandian and Laurentide ice sheets alters the mid-latitude circulation pattern and shifts the 850 hPa jet stream southward, leading to decreased precipitation over the ice-covered region but increased precipitation over its southern periphery.
Similar to the PMIP3 models, the NESM v1 simulates higher hydrological sensitivity during the LGM period than the hydrological sensitivity in future projection. Results of the single forcing experiments suggest that the higher hydrological sensitivity could be primarily attributed to the LSC forcing in the global mean sense, because it amplifies the sensitivity in the convergence zones of both hemispheres where precipitation and hydrological sensitivity are large. Detailed analysis reveals that the dynamic effect related to circulation change is the key.