PMIP4 experiments using MIROC-ES2L Earth System Model

. Following the protocol of the fourth phase of the Paleoclimate Modeling Intercomparison Project (PMIP4), we performed numerical experiments targeting distinctive past time periods using the Model for Interdisciplinary Research on Climate, Earth System version 2 for Long-term simulations (MIROC-ES2L), which is an Earth System Model. Setup and basic performance of the experiments are presented. The Last Glacial Maximum was one of the most extreme climate states during the Quaternary and conducting numerical 15 modeling experiments of this period has long been a challenge for the paleoclimate community. We conducted a Last Glacial Maximum experiment with a long spin-up of nearly 9,000 years. Globally, there was reasonable agreement between the anomalies relative to present day derived from model climatology and those derived from proxy data archives while some regional discrepancies remained.


Introduction
Using climate models to simulate past climate provides unique opportunities to evaluate models' projections of future climate.
The Paleoclimate Modeling Intercomparison Project (PMIP) began in the early 1990s (Joussaume et al., 1999).Since then, the paleoclimate community has continued to expand their research to include more time periods and events.With the increase of computational power, models of higher complexity are used to make future projections (Kawamiya et al., 2020).Phase 3 of PMIP was endorsed by phase 5 of the Coupled Model Intercomparison Project (CMIP5; Braconnot et al., 2012), and PMIP is now in its fourth phase (PMIP4; Kageyama et al., 2018).The proposed PMIP4 experiments cover a wide range of time periods, including the Last Glacial Maximum (LGM; 21,000 years before present), two interglacials (6,000 and 127,000 years before present), the last millennium (LM), the mid-Pliocene, and many non-CMIP time periods.
The Quaternary is characterized by cyclic climate change with long glacials and short interglacials that have been recorded in various paleo-proxy records such as ice cores (Jouzel et al., 2007;Dome Fuji Ice Core Project Members, 2017), ocean sediment cores (Weldeab et al., 2007), loess records (Maher et al., 2010), and terrestrial fossils (Bartlein et al., 2011).The LGM refers to the period when global ice volume reached its maximum.It was also one of the coldest periods of the Quaternary.
Since the beginning of PMIP, attention has been drawn toward the LGM, which was one of the extreme periods in the glacialinterglacial cycles of the Quaternary (Joussaume and Taylor, 2000;Braconnot et al., 2007;Kageyama et al., 2020), and also the most recent period during which global coverage of the continental ice sheets was at its maximum and greenhouse gas (GHG) levels were at a minimum.
As LGM cooling relative to the preindustrial (PI) experiment over the tropics is at a comparable level to equilibrium climate sensitivity (ECS), LGM modeling can provide useful information to constrain climate sensitivity for projections of future climate (Annan et al., 2006;Renoult et al., 2020).Intercomparison studies of proxy-based reconstructions of climate variables and model output continue to be conducted (Braconnot et al., 2007;Bartlein et al., 2011;Kageyama et al., 2020).They report good agreement between model output and proxy data for temperature and sea surface temperature (SST) anomalies over the low latitudes (Otto-Bliesner et al., 2009;Hargreaves et al., 2013); however, the tendency for models to underestimate cooling over Greenland remains (Masson-Delmotte et al., 2006).Models have difficulty in reproducing the weakened Atlantic meridional overturning circulation (AMOC) of the LGM (Weber, 2007;Brady et al. 2013;Muglia and Schmittner, 2015;Marzocchi and Jansen, 2017), which might influence the underestimation of cooling.The dust deposition was several to tens of times higher at LGM (Lambert et al. 2008, Lamy et al. 2014, Dome Fuji Ice Core Project members 2017), but was difficult to reproduce by LGM experiments; to reproduce the dust abundance at LGM, we need to assume glaciogenic dust (Mahowald et al. 2006, Ohgaito et al. 2018), or assuming an erodibility map (Albani et al. 2014).And an erodibility map was formally introduced in PMIP4 ( Kageyama et al. 2017), in addition to the dust emission that is simulated in non-Paleo simulations.In Ohgaito 2018, they showed that sufficient dust loading affects the temperature around Antarctica.
The interglacial periods of 6,000 and 127,000 years before present were characterized by differences in solar radiation at the top of the atmosphere caused by orbital states that were different from those of the present day (Brierley et al., 2020;Otto-Bliesner et al., 2020), resulting in seasonalities that were different from the PI period (1850 CE).As it was in the recent past and because various paleo-proxy records are available (Ritchie et al., 1985;Drake et al., 2011;Hely et al., 2014;Tierney et al. 2013Tierney et al. , 2017)), the interglacial period of 6,000 years before present was the only interglacial included in earlier phases of PMIP (Braconnot et al., 2007;Ohgaito andAbe-Ouchi, 2007, 2009;Ohgaito et al., 2013).Following efforts to collect paleo-proxy data (Otto-Bliesner et al., 2001;Lunt et al., 2013;Capron et al., 2014Capron et al., , 2017)), it is now also possible to conduct the same experiment for 127,000 years before present; experiments on this interglacial period have the advantage of strong seasonality in the Northern Hemisphere (NH).The insolation anomaly at 127,000 years before present was larger than that at 6,000 years before present, and the stronger insolation at 127,000 years before present during boreal summer modulated the temperature and circulation for that time period (Lunt et al., 2013).The role of vegetation coupling has been discussed intensively as studies report that vegetation enhances warming in the NH (O'ishi and Abe-Ouchi, 2011).However, models have been unable to reproduce the quantitative changes recorded in proxy data (McKay et al., 2011, Capron et al., 2014, Hoffman et al., 2017).
As the LM is the most recent period prior to the PI period, there are vast amounts and varieties of data available from exact times in proxy records (PAGES2k-PMIP3 group, 2015;Luterbacher et al., 2016;Gagen et al., 2016) and in the literature (Pfister and Brazdil, 2006;Xoplaki et al., 2016;Camenisch et al., 2016).In earlier numerical paleoclimate studies, simple models were used to conduct transient experiments over periods of 1,000 years (Crowley, 2000;Goosse et al., 2005).With the increase of computational power, simulations using coupled Atmosphere-Ocean General Circulation Models (AOGCMs) and/or comprehensive Earth System Models (ESMs) became standard (Kawamiya et al., 2020).Coordination of LM experiments began under PMIP3 (Schmidt et al., 2011), and multiple AOGCMs and ESMs have been used to perform LM experiments.One of the important questions for LM experiments is whether climate variabilities stem from internal variability or forced responses.Atwood et al. (2016) decomposed the forcing of the LM experiment and concluded that cooling during the Little Ice Age (LIA; 1450-1850 CE) was largely driven by volcanic eruptions.PAGES2k (2015) summarized reconstruction-model intercomparisons and reported that the agreement between model and proxy-based reconstructions is better in the high latitudes in the NH and worse in the Southern Hemisphere.Historical (HIST;1850-2014 CE) CE) and LM experiments are intrinsically different from the other PMIP4 time-slice experiments discussed in this paper.They are timevarying experiments that follow the same method used in the historical experiment in CMIP6.Hence, the LM experiment is closely aligned with the scientific focus of other endorsed MIPs, such as comparison of climatic response to volcanic forcing (VolMIP; Zanchettin et al., 2016) and land use (LUMIP; Lawrence et al., 2016).
Using the Model for Interdisciplinary Research on Climate, Earth System version 2 for Long-term simulations (MIROC-ES2L), we performed numerical experiments targeting distinctive time periods.These simulations included the LGM, the 6ka and 127ka, and the LM experiments.The MIROC-ES2L is an ESM that contains atmosphere, ocean, land, and ocean and land biogeochemical cycles (Hajima et al. 2020), which has been developed recently to contribute to CMIP6 and the United Nations Intergovernmental Panel on Climate Change Sixth Assessment Report.
The model is presented in Sect. 2 and the experimental setup and spin-up procedures are explained in Sect.3. Basic climate states from the experiments are presented in Sect.4, and conclusions and an outlook are discussed in Sect. 5.

Model
The MIROC-ES2L is an ESM developed for CMIP6 (Hajima et al., 2020) and its physical core comprises atmosphere, ocean, and land modules; variables are exchanged via a flux coupler (Fig. 1).The AOGCM components are the same as those in Tatebe et al. (2018).The physical ocean and land modules are coupled with the land ecosystem model VISIT-e (Ito and Inatomi, 2012) and the ocean biogeochemical model OECO2 with a nutrient-phytoplankton-zooplankton-detritus type representation of the ecosystem.The ecosystem modules can simulate global carbon and nitrogen cycles explicitly.As the carbon and nitrogen in the atmosphere are prescribed to predetermined values in each experimental setting in this study, carbon and nitrogen variables calculated in OECO2 and VISIT-e are not returned to the atmosphere.Distribution of plant functional types (PFTs) is prescribed because VISIT-e is not a dynamic vegetation model.Dynamics of aerosols are calculated by an online aerosol module, SPRINTARS (Takemura et al., 2000(Takemura et al., , 2005(Takemura et al., , 2009)).
Horizontal resolution of the atmosphere is set to T42 spectral truncation.Vertical resolution is 40 levels up to 3 hPa.The ocean component has tripolar horizontal coordinates, with two poles in the NH that are located over land to avoid singularity over ocean grids.Horizontal resolution of the ocean is 1° longitude and varies from 0.5° latitude around the Equator to 1° latitude over the mid-latitudes.Vertical resolution of the ocean is 62 layers with a hybrid sigma-z coordinate.
Using this model, various types of CMIP6 experiments have been performed.These include all of the Diagnostic, Evaluation and Characterization of Klima (DECK) experiments, the historical experiment of CMIP6 (from 1850 CE to 2014 CE), and the endorsed MIP experiments.The ECS of this model version is 2.66 K by Gregory's method (Gregory et al. 2004).The identical model version was used for all the experiments in this study.

Setup and spin-up of preindustrial (PI) control experiment
The PI control experiment is the reference experiment of all the paleoclimate experiments.It is identical to the piControl experiment in CMIP6 (Eyring et al., 2016) and the experimental configuration of PI in MIROC-ES2L is described in detail in Hajima et al. (2020).Levels of GHGs were set following the protocol of CMIP6: CO2, CH4, and N2O were set to 284.725 ppm, 808.25 ppb, and 273.02 ppb, respectively (Table 1).The PFTs in PI are inherited from MIROC-ESM (Watanabe et al. 2011), which was based on Ramankutty and Foley (1999) (Fig. 2(c)).A description of each PFT is given in the caption of Fig. 2(c).
Topography is defined from GTOPO30 (Fig. 2(e)).The experiment was run for more than 9,000 model years during the course of model development and the final drift of the global mean surface air temperature was −4.79 × 10 −5 °C yr −1 for the final 500 years (Hajima et al. 2020).Model output from this period was submitted to CMIP6 and the climatology of this period is used for the analyses in this study.

Setup and spin-up of Last Glacial Maximum (LGM) experiment
We performed the LGM experiment following PMIP4 protocol (Kageyama et al., 2017).A long spin-up is essential because of the considerable differences between LGM and present-day conditions.Hence, before model development was finalized, we started spinning up using the physical core (AOGCM) of MIROC-ES2L (Tatebe et al., 2018).Spin-up started by reducing CO2 (Bereiter et al., 2015), CH4 (Loulergue et al., 2008), and N2O (Schilt et al., 2010) levels from PI to LGM values (Table 1).Global mean air temperature gradually reached quasi-equilibrium.After integration for 2,640 model years, the land-sea mask, ice sheets, altitude (from ICE-6G_C, as presented in Peltier et al. (2015) (Fig. 2(b, d, and f))), river courses, and Earth's orbit (Berger, 1978) (Fig. 3(a)) were changed from PI to LGM conditions step by step.The total spin-up time was 6,760 model years (Figs. 4 and 5).The LGM PFTs were created based on the PI PFTs with the ice sheet grids defined by ICE-6G_C, and nearby PFTs were diverted to non-ice sheet land (exposed continental shelves) expanded from PI.In this process, the biggest technical difficulty occurred when the LGM land-sea mask and ice sheet albedo were introduced ((2) of Fig. 5).As mentioned in Appendix, the LGM spin-up failed at year 3305 and the simulation could not proceed because a grid reached the upper limit of the model's sea ice thickness on the north coast of North America.Figure A1 shows the time evolution of winter sea ice thickness along the north coast of North America during the spin-up.In the present-day condition, the Beaufort gyre prevails on the north coast of North America, with southward winds pushing sea ice toward land.In addition to this, the sea ice thickness has increased too much due to the cold conditions.When we changed the elevation to the LGM condition, which was left at PI ((3) of Fig. 5), the circulation field changed due to the appearance of the high altitude Laurentide ice sheet, and the north coast of North America had neutral north-south winds on average, which would eliminate the sea ice build-up phenomenon and allowed the experiment to continue.
As the development of MIROC-ES2L was finalized during LGM spin-up, conditions in the 6,760 th model year of the spin-up were used to initiate the LGM experiment in MIROC-ES2L.In this conversion procedure, the offline terrestrial module was spun up for 40,000 model years until quasi-stability was reached and the end state was used in the LGM experiment in MIROC-ES2L.This was followed by spinning up the main MIROC-ES2L experiment for a further 100 years.Ocean salinity (1 Practical Salinity Unit (PSU) was added globally) and an erodibility map (addressing dust emission under LGM conditions, as proposed by Albani et al. (2014Albani et al. ( , 2016)) were introduced.The erodibility map specifies low latitudes as deserts and mid-to high latitudes as tundra (Fig. 2(d)).Land and ocean ecosystem models were spun up offline for 40,000 and 3,000 model years, respectively, on the basis of the physical conditions created by MIROC-ES2L.Land and ocean biogeochemical states at the end of the offline spin-up were used to initialize the LGM experiment in MIROC-ES2L.The LGM experiment was run for a further 1,800 years until it eventually reached quasi-equilibrium.Surface air temperature of the final 500 model years showed a trend of 0.0002 °C yr −1 .Model output from the final 100 years was submitted to PMIP4-CMIP6.

Setup and spin-up of the two interglacial experiments
The 6ka and 127ka experiments were spun up following the protocol outlined in Otto-Bliesner et al. (2017).The specified GHGs and orbital parameters were as listed in Table 1.The main difference between these periods and the PI period was the change in insolation attributable to Earth's orbit, as shown in Fig. 3(b and c), where seasonality was amplified in the NH and diminished in the Southern Hemisphere.
Starting from PI, the 6ka experiment was integrated for 1,500 model years and the 127ka experiment was integrated for 1,550 model years (Fig. 4(b and c)).After the long spin-up, the final 100 years of the simulations were selected as the formal products to be submitted to CMIP6 and for analyses in this study.The 127ka experiment is identical to the LIG experiment in O'ishi et al. (2020).

Setup and spin-up of the Last Millennium and historical experiments
We performed the LM experiment following the protocol of Jungclaus et al. (2017).The experiment was forced with timevarying total solar irradiance (Shapiro et al., 2011;Vieira et al., 2011;Wu et al., 2017), orbit (Berger, 1978), GHGs (Meinshousen et al., 2017), volcanic eruptions (Sigl, 2015;Tooney and Sigl, 2017), ozone, and land use change (Hurtte et al., 2016) (Table 1).The experiment followed the same procedure as that of the historical experiment in CMIP6 (Eyring et al., 2016).From PI, the model was run under constant forcing from 850 CE for 200 model years (Fig. 4 (d)).The end state of the spin-up was used to initialize the time-varying LM experiment, which was conducted from 850 CE to 1850 CE.We performed a HIST experiment following CMIP6 protocol (Eyring et al., 2016).The end state of the LM experiment was used to initialize HIST, which was run until 2014 CE.Hajima et al. (2019) analyzed basic model performance for the present and indicated that MIROC-ES2L is a state-of-the-art ESM that is able to reproduce mean climatology reasonably well.Global annual mean air temperature at 2 m height is 14.99 °C and the peak value of the annual mean AMOC is 15.3 Sv, which falls within the range of reasonable estimates (Frajka-Williams et al., 2019).The peak AMOC is defined as the peak value of the area from 15°-60°N and from 900-3300 m depth.Model SST presents a reasonable global distribution but has positive bias over the Southern Ocean, which leads to underestimation of the extent of Antarctic sea ice.

LGM mean climate
Relative to PI, the final 100 years of the LGM has a global mean surface air temperature anomaly of −4.4 °C (Figs.5(a) and 6(a)) and a tropical air temperature anomaly of approximately −2 °C, consistent with values derived from paleo-proxy archives (MARGO project members, 2009;Bartlein et al., 2011).Borehole thermometry suggested a temperature anomaly of the LGM relative to PI over Eastern Antarctica of −7 to −10 °C (Stenni et al., 2010;Uemura et al., 2012).The temperature anomaly in the model is approximately −6.0 °C, suggesting that cooling in the model is weak.For central Greenland, borehole thermometry suggested a temperature anomaly of −21 to −25 °C (Cuffey et al., 1995;Johnsen et al., 1995;Dahl-Jensen et al., 1998), whereas the model temperature anomaly is −11.1 °C.The large discrepancy between ice core data and model output could partly be attributed to issues related to the modeling of the thermohaline state of the ocean (McManus et al., 2004;Curry and Oppo, 2005).It is well known that numerical models have difficulty in reproducing the sluggish thermohaline circulation of the LGM (Otto-Bliesner et al., 2007;Muglia and Schmittner, 2015;Sherriff-Tadano et al., 2017) that is suggested in proxy data (Lynch-Stieglitz et al., 2007;Hesse et al., 2011).In our experiment, the peak value of the annual mean AMOC of the LGM is 21.0 Sv (Figs. 5(b) and 7), which is higher than that of PI.To address this issue, we will continue the experiment to identify the components that contribute to global cooling and those that contribute to cooling over the polar regions.The sharp gradient shown by the contour lines around 32°N would be caused by a strong and deep westerly boundary current and associated strong upwelling (Brady et al., 2013), which can be seen in the previous LGM modelling studies having strong AMOC (Brady et al. 2013, Muglia and Schmittner 2015, Sherriff-Tadano et al. 2018).
Figure 6(b) shows the net precipitation anomaly relative to PI.Total precipitation is 1063 mm yr −1 for LGM and 1166 mm yr −1 for PI.Consistent with Bartlein et al. (2011), model precipitation has a general tendency to be lower for the LGM than for the PI period because the lower SSTs and colder climate of the LGM result in a weaker hydrological cycle, which is also shown in the weakened zonal water vapor transport (Fig. S1(b)).Large reductions in precipitation relative to PI are found in areas that were covered by ice sheets during the LGM but were no longer ice-covered at PI, i.e., the areas covered by the Laurentide and Fennoscandian ice sheets.These large anomalies would be associated with the higher altitude of the ice surface relative to the ground surface when the ice sheets have disappeared.In the northern North Atlantic Ocean, large anomalies are associated with the southward expansion of sea ice during the LGM.
Anomalies of zonal mean oceanic potential temperature and salinity are shown in Fig. S2(a and d).In the Southern Ocean, proxy data (Adkins et al., 2002) suggested anomalies of −2 °C and +2.5 PSU at around 3600 m depth.The underestimation might be attributed to too little sea ice formation, which would be related to a warm bias of the Southern Ocean (Hajima et al., 2020).In contrast, in the North Atlantic Ocean, the anomaly of salinity agrees with the proxy data (Adkins et al., 2002), whereas the temperature anomaly (−1 to −2 °C) is underestimated (−4 to −5 °C; Adkins et al. (2002)).A temperature that is too warm could possibly be attributed to a high state of AMOC.
In Fig. 8, we compare the dust deposition fluxes with data archives (Kohfeld et al., 2013, Albani et al., 2014).The distribution is also shown in Fig. S3.The model shows general consistency with the data archives globally, with positive bias in the PI over Antarctica and Greenland, and values that are insufficiently high in the Gobi and Taklamakan regions.The LGM shows better representation of the proxy data than the PI, with reasonable fluxes over Antarctica.However, it underestimates the high values that are abundant in the East Asian region and the high dust fluxes in North America.The LGM dust fluxes are shown in Fig. S3(c) as a ratio of the LGM dust fluxes to those of the PI.The ratio is generally well represented globally, but the ratio is underestimated in South America and in regions of the South Atlantic downstream of the wind.The reason for the underestimation, as mentioned above, is probably overestimation of South American dust emission in PI.
Figure 9 shows the export production anomaly of the oceanic ecosystem of the LGM relative to PI with paleo-proxy data (Kohfeld et al., 2013) superimposed on model output.As proxy data provide qualitative information rather than quantitative assessments, comparisons between model output and proxy data can only be used to evaluate the accuracy of the general direction of the model anomaly.Positive model anomalies over the low-mid-latitudes are consistent with the proxy data.
Negative anomalies over the high latitudes can be understood to reflect sea ice expansion during the LGM.Sea ice expansion during the LGM is underestimated in the model (Crosta et al., 1998) and could result in negative anomalies around Antarctica that have smaller absolute values than those indicated by proxy data (Fig. 9).
An excessively weak positive anomaly around 40°-50°S in LGM could be the result of dust emission being too high over South America in the PI experiment (Fig. 8).Mean global terrestrial gross primary production in the LGM experiment is 65% of that in the PI experiment, which is consistent with the estimates of Prentice et al. (2011) 2 and Fig. S4).The simulated glacial ocean is therefore unable to explain the glacial-interglacial drawdown of atmospheric CO2, which is similar to previous modeling studies (Buchanan et al., 2016).It should be noted that the effects of burial-nutrient feedback and carbonate compensation on the oceanic carbon cycle are not considered in this simulation because MIROC-ES2L does not include a sediment module.
As the dissolved oxygen cycle is the mirror image of the biological carbon cycle, reconstructed oxygen change is useful to constrain the respired carbon accumulation.We compared modeled oxygen changes from PI to LGM with qualitative and quantitative proxy data (Jaccard et al., 2016;Durand et al., 2018;Schmiedl and Mackensen, 2006;Hoogakker et al., 2015Hoogakker et al., , 2018;;Gottschalk et al., 2016;Lu et al., 2016;Bunzel et al., 2017;Umling and Thunell, 2018).The combination of cooler SST and enhanced AMOC increases the oxygen content by approximately 10 mmol m −3 in both the upper and the deep ocean (Table 2).The simulated increases in oxygen are in reasonable agreement with the proxy data for the upper ocean, but contrast the proxy data for the deep ocean, which show a decrease in oxygen of more than 30 mmol m −3 (Fig. S4).The model-proxy disagreements of deep ocean oxygen change result in underestimation of the accumulated respired carbon.
The simulated deep-water-mass age is younger during the LGM than during the PI by approximately 60 years (Table 2), indicating an increase in ventilation due to enhanced AMOC.However, proxy data show an increase in water-mass age of more than 1,000 years, suggesting reduced ventilation and weaker AMOC (Burke and Robinson, 2012;Curry and Oppo, 2005).
Enhanced ventilation supplies oxygen-rich surface water to the deep ocean and simultaneously releases carbon accumulated in the deep water to the atmosphere.Therefore, we attribute the model-proxy disagreements of deep ocean oxygen change and underestimation of respired carbon accumulation to overestimation of ventilation.Our results suggest that weaker AMOC is required for reproducing the respired carbon accumulation and deoxygenation in the glacial deep ocean.

Mean changes in interglacial experiments
Figures 10 and 11 show air temperature and precipitation anomalies of 6ka and 127ka relative to PI.Because of the marked changes in seasonality in these time periods, adjustment of the calendar was applied and average anomalies were calculated for June-August (JJA) and December-February (DJF).Air temperature anomalies are positive over continental interiors in the mid-high latitudes in JJA for both 6ka (Fig. 10(a)) and 127ka (Fig. 11(a)) because of changes in shortwave radiation forcing.
Compared with 6ka, stronger shortwave forcing results in larger air temperature anomalies in 127ka (Fig. 3).
Precipitation anomalies suggest that, relative to PI, boreal summer monsoons were stronger (Figs. 10(b) and 11(b)) and austral summer monsoons were weaker during the interglacials (Figs. 10(d) and 11(d)).The precipitation anomaly over the Sahel region suggests that the reduction of desert area during the interglacials relative to PI is smaller in the model than that suggested by proxy data (Petit-Maire, 1999;Castaneda, 2009;Tierney et al., 2017;Drake et al., 2011;Hely et al., 2014), which is a mismatch that has been persistent through many modeling efforts (Braconnot et al., 2001(Braconnot et al., , 2007)).The zonal water vapor transport is shown in Fig. S1(c and d).Relative to PI, more water vapor is transported to North Africa, and the amplitude is more pronounced in 127ka than 6ka following the radiation anomaly.
Variations of temperature and precipitation anomalies with season and latitude are shown in four Hovmöller diagrams in Fig.
12. The temperature anomaly responds to changes in solar radiation with a lag of approximately one month (Figs.3(b and c) and 12(a and b)), which could be a consequence of the slow thermal response of the ocean surface.Precipitation anomalies of 6ka and 127ka, relative to PI, exhibit a northward shift and enhancement during boreal summer in the NH (Fig. 12(c and d)) which is consistent with Figs.10(b) and 11.Fig. S2 shows the anomaly of the zonal mean oceanic potential temperature and salinity of 6ka and 127ka relative to PI.
Surface cooling is consistent with Fig. 12 and freshening would be result of a more active water circulation in the NH.Strong freshening around 32°N in Fig. S2(f) is attributed to low salinity in the Mediterranean Sea.The LIA is reasonably well expressed in the NH mean, but the warming during the Medieval Climate Anomaly (MCA; 950-1250 CE), which is suggested by proxy data, is underestimated by the model.The difference between the NH mean temperature of the LIA and that of the MCA is −0.1 °C, which is not statistically significant in the Student's t-test.

Last Millennium and historical transient variabilities
The HIST experiment was run for the period between 1850 and 2014 CE.

Outlook and conclusions
Using MIROC-ES2L, an ESM that has recently been developed for CMIP6, we performed numerical experiments to examine the paleoclimate during several time periods and one historical experiment that was initiated from the LM experiment.
Globally, there was reasonable agreement between the climate states described by the model and those derived from proxy data, while some regional discrepancies remained.In this section, we summarize the results and explore the possible causes of the discrepancies.
From PI, LGM conditions were introduced step by step into MIROC-ES2L and the LGM spin-up experiment was run for approximately 9,000 model years.The temperature anomaly of LGM relative to PI over the tropics is negative and there is general quantitative agreement between the anomaly derived from the model and that from proxy data (Bartlein et al., 2011;MARGO project members, 2009).This could be useful in constraining future projections, given that Annan et al. (2005), and Renoult et al. (2020) revealed the correlation between tropical cooling at LGM and ECS.It has been pointed out in the United Nations Intergovernmental Panel on Climate Change (2013) Fifth Assessment Report that the cooling over Greenland during the LGM relative to PI is underestimated in the models.This could be attributed to strong AMOC in the models, which leads to an estimate of sea ice expansion over the northern Atlantic Ocean that is lower than suggested by proxy data.The anticorrelation between sea ice expansion and AMOC is known from observations (Boehm et al. 2015) and modeling (Peltier and Vettretti, 2014).Intrusion of Antarctic bottom water into the Atlantic Basin is very weak in MIROC-ES2L, even in the PI experiment (Tatebe et al., 2018).Insufficient abyssal flow into the Atlantic Basin could be partly caused by the low resolution of the ocean component.Detailed analyses on the representation of atmospheric circulations would be necessary for further investigation.Model representation of the Southern Ocean might influence the distribution of CO2 between the atmosphere and the ocean (Moore et al., 2000).Anomalies associated with topography might be obscured by the low horizontal resolution of the model, resulting in discrepancies between climate states in the model and those derived from proxy data.Cooling of Eastern Antarctica during the LGM relative to PI, which is suggested by ice core data (−7 to −10 °C), is underestimated by this model (−6 °C), as explained in Sect.4.2.This could be partly attributed to the positive SST bias over the Southern Ocean in the model at PI and subsequent underestimation of sea ice expansion.PMIP model analyses (Otto-Bliesner et al. 2007, Marozzochi andJansen 2017) also suggested the correlation of AMOC and sea ice coverage.
There is reasonable agreement between dust flux from the LGM experiment and that suggested by proxy data.However, the PI experiment overestimates dust emissions from South America.Thus, the change in dust emission between LGM and PI is likely to be underestimated in the model, leading to underestimates of LGM anomalies relative to PI for climate (Ohgaito et al., 2018) and ecosystem activity in the Southern Ocean (Yamamoto et al., 2019).Further studies will be necessary to investigate the impact of representing reasonable dust emission and loading on climate.
We prescribed conventional land PFTs in the LGM experiment, as discussed in Sect.3.2.In models that comprise a coupled dynamic vegetation model, climate states would be altered through biophysical feedback because of changes in vegetation cover (O'ishi and Abe-Ouchi, 2013).
The LGM experiment was also performed in the previous phase of PMIP using MIROC-ESM (Sueyoshi et al., 2013), which is the previous version of MIROC-ES2L.Owing to differences in forcing (mainly in terms of GHGs and ice sheets) and spinup procedures, we are unable to make direct comparisons of the experiments conducted using the two versions of the model.However, there is a general tendency of the PMIP4 model to simulate less cooling at LGM relative to PI, which is a tendency that was also identified by Kageyama et al. (2020) in their comparison of LGM experiments from different versions of PMIP.
Further sensitivity experiments using different boundary conditions could be helpful for identifying causes of this discrepancy.
Although we conducted long spin-up for the LGM experiment, the abyssal salinity and oceanic temperature are not representative of the structure suggested by proxy data.This discrepancy might reflect model biases, e.g., SST bias, and/or difficulty in representing the AMOC state of the LGM.
The two interglacial experiments (6ka and 127ka) include different orbital parameters and GHG levels, and had long spin-up times that exceeded 1,400 years.Results showed warming over NH continents during boreal summer relative to PI, consistent with the direction of change suggested by proxy archives (Bartlein et al., 2011;Turney and Jones, 2010); however, the model underestimates the amount of warming.The discrepancy could be reduced by improving the experimental setup, such as replacing the prescription of PFTs by a process that could produce PFTs that are closer to the real conditions of the periods.
The vegetation coupling greatly improves the representation of the warmings shown by proxies at the Arctic Ocean margin (O'ishi and Abe-Ouchi, 2011, O'ishi et al. in press CP).On the other hand, some inconsistency remains in inland areas such as inner Eurasia.Pfeiffer and Lohmann (2016) suggested that we need to take into account the uncertainty of the times of the proxy data.
Compared with PI, temperature over the tropics is generally lower in the 6ka experiment, which may partly be inconsistent with proxy data archives (Bartlein et al., 2011, Kaufman et al., 2020).Because the uncertainty of the marine proxies is so large (Hessler et al., 2014) that it is unlikely to be possible to identify positive or negative SST anomalies at 6ka, we exclude them from our assessment here.However, cooling over the tropics could be considered a reasonable and direct response to net negative solar forcing.Thus, the discrepancy between the model and proxy data suggests that feedbacks that might play a role in the modeling of climate change are missing in the current model.Further improvement and expansion of the model would be necessary.
The precipitation anomaly shows a northward shift of peak precipitation in boreal summer in the NH.The precipitation anomaly over the Sahara Desert in the model is still smaller than that suggested by the proxy data archive, which is a mismatch that has been persistent throughout the long history of PMIP.It might be necessary to include new processes to maintain higher soil moisture in the interglacials (Hopcroft et al., 2017).
The LM experiment performed in this study showed clear responses of global temperature to huge volcanic eruptions, while the responses of global temperature to other forcings were unclear.Responses to external forcings except volcanos are likely to be small in comparison with internal variabilities.
The difference between model NH mean temperature at the LIA and that at the MCA is too small to be statistically significant.
However, earlier studies suggested that signals might be more pronounced at regional scales (PAGES2k, 2015;Fernandez-Donado et al., 2013); thus, further investigations regarding regional scales would be necessary.The HIST experiment was initiated from the end of the LM experiment and it produced time series of global temperatures that are similar to those from the other historical experiments that were initiated from PI (Hajima et al., 2020).This suggests that the initial conditions used for the standard historical experiment in CMIP6 are appropriate for the simulation of global temperatures in the industrial era.
Sensitivity experiments using different boundary conditions would be useful for identifying causalities to obtain more details in future analyses.
obtained using a dynamic global vegetation model.The dissolved inorganic carbon content of the ocean in LGM is 629 Pg C less than in PI.Lowering the atmospheric CO2 to 190 ppm and strengthening of the overturning circulation lead to extraction of a large amount of carbon out of the ocean.Conversely, increased solubility owing to cooling and enhanced biological carbon export because of increases in nutrient and iron supply from the ocean interior and dust lead to accumulation of carbon within the ocean.The former effects mainly contribute to carbon redistribution, resulting in reduction in the carbon content of 239 Pg C in the upper ocean and 390 Pg C in the deep ocean (Table

Figure 13
Figure 13(a) shows the time series of annual mean NH air temperature of LM and HIST.Sharp cooling events are clear responses to huge volcanic eruptions.The effect of solar forcing on annual mean temperature is unclear, probably because the signals are small in comparison with internal variability.
Figure 13(b) shows the output of 30 ensembles of MIROC-ES2L historical experiments submitted to CMIP6 and data from HadCRUT4 (Morice et al., 2012), which is scaled at the mean value of 1960-1989.Centennial variabilities of the NH mean temperature in HIST and in the CMIP6 historical experiments consistently show a positive trend during the first half of the 20 th century, followed by a cooling trend until 1970 and then subsequent warming.In comparison with the standard historical experiments, HIST has a less positive bias (Fig. 13(c)).

Figure 5 :
Figure 5: Time series for spin-up and submitted period (final 100 years) of LGM experiment and PI as a reference for (a) global mean air temperature at 2 m height and (b) peak values of annual mean AMOC.Gray line (0) denotes PI value.Blue line (1): the experiment with only GHG levels set to LGM values.Light blue line (2): with the land-sea distribution and land PFTs hanged to 835

Figure 8 :Figure 9 :
Figure 8: Dust deposition from model output and derived from proxy archives (Kohfeld et al., 2013; Albani et al., 2014) for (a) PI and (b) LGM (g m −2 yr −1 ).Colors represent the locations of the proxy data, as explained in the legend in the figure.Crosses, circles, and diamonds represent terrestrial, marine, and ice core data.

Figure 12 :
Figure 12: Hovmöller diagrams for (a) air temperature anomaly at 2 m height (°C) for 6ka relative to PI, (b) air temperature anomaly at 2 m height (°C) for 127ka relative to PI, (c) precipitation anomaly (mm d −1 ) for 6ka relative to PI, and (d) precipitation anomaly (mm d −1 ) for 127ka relative to PI. Calendar adjustments are applied.

Figure 13 :
Figure 13: (a) Annual mean air temperature (°C) averaged over the Northern Hemisphere from the LM and HIST experiments (red: annual mean, dark red: 10 year running mean) and observational data from HadCRUT4 (black).(b) Annual mean air temperature (°C) averaged over the Northern Hemisphere from 1850-2014 for HIST and 30 ensemble members of the historical experiments (gray) and HadCRUT4 (black).HadCRUT4 is scaled for the period from 1961-1990 (period shaded light green).(c) Histogram of deviations from HadCRUT4 shown in (b), averaged for every five-year mean during 1850-1949, counted in 0.1 °C 870

Figure A1 :
Figure A1: Average sea ice thickness (m) from January to March for the first half of the spin-up of the LGM experiment, averaged 150° -180° E, 70° -75° N. The blue line: step 1, light blue: step 2, and yellow-green: step 3 in Figure 5, respectively.For simplicity,