PMIP4/CMIP6 Last Interglacial simulations using different versions of MIROC, with and without vegetation feedback

We carry out a Last Interglacial (LIG) experiment, named lig127k and a Tier1 experiment of PMIP4/CMIP6, using three versions of the MIROC model, MIROC4m, MIROC4m-LPJ and MIROC-ES2L. The results are compared with reconstructions from climate proxy data. All models show summer warming over northern high latitude land, reflecting the 10 differences between the distributions of the LIG and present-day solar irradiance. Only MIROC4m-LPJ, which includes dynamical vegetation feedback from the change in vegetation distribution, shows warming signals, even for the annual mean, at northern high latitudes, as shown by proxy data. However, the latest Earth System Model (ESM) of MIROC, MIROC-ES2L, in which there is only a partial vegetation effect through the leaf area index, shows no change or even annual cooling over large parts of the northern hemisphere. Results from the series of experiments show that the inclusion of vegetation feedback 15 is necessary for the reproduction of the strong annual warming over land at northern high latitudes. The LIG experimental results show that the warming predicted by models is still underestimated, even with dynamical vegetation, compared to reconstructions from proxy data, suggesting that further investigation and improvement to the climate feedback mechanism are needed.


Introduction 20
The Last Interglacial (LIG, 130ka-116ka) is referred to as the warmest period in the recent glacial-interglacial cycle (NGRIP members 2004, Overpeck et al. 2006. The most important characteristic of the LIG is the strong summer solar irradiance in the northern hemisphere due to the difference in the earth's orbit during that period and that of the present-day (Berger 1978).
For example, 127,000 years ago, summer solar irradiance was more than 50W/m 2 larger compared to that of the present-day ( Figure 1). Paleo-evidence shows strong warming at northern high latitudes in response to this different spatial and temporal 25 pattern of solar irradiance (Otto-Bliesner et al. 2006, Lund et al. 2013. Northward shift of boreal treeline due to warming is also indicated by proxies (LIGA members 1991, Edwards et al. 2003. Sea level rise due to warming has also been pointed out, with contributions from the mass balance change in the Greenland and the Antarctic ice sheets (Robinson et al. 2011, Born and Nisancioglu, 2012, Quiquet et al. 2013, Stone et al. 2013. https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License.
The Paleoclimate Modelling Intercomparison Project (PMIP) coordinates the cooperation and comparison between modelling 30 and data of the past (Braconnot et al. 2000, Braconnot et al. 2007, Braconnot et al. 2012. The LIG is one of targeted periods as well as the mid-Holocene and the Last Glacial Maximum . For the simulation of the past periods, PMIP provides protocols with common settings which should be applied to participating models. In the present study, we apply the LIG boundary conditions provided by the PMIP4 to three different version of atmosphere-ocean coupled general circulation models (AOGCMs) which belong to the Model for Interdisciplinary Research on Climate (MIROC) family and 35 compare results with the pre-industrial simulations, focusing on the different treatment of vegetation among the three models.

MIROC4m
An AOGCM, MIROC4m, is based on MIROC3.2 which contributed to the fourth assessment report (AR4; Meehl et al. 2007 Hasumi 2006). The AGCM solves the primitive equations on a sphere using a spectral method (Numaguti et al. 1997).
The model resolution of atmosphere component is T42 with vertical 20 layers. The land submodel is the Minimal Advanced Treatments of Surface Interaction and Runoff (MATSIRO; Takata et al. 2003). Vegetation is prescribed to a satellite based 45 present-day distribution. The OGCM solves the primitive equation on a sphere, where the Boussinesq and hydrostatic approximations are adopted (Hasumi 2006). The horizontal resolution is ~1.4° in longitude and 0.56°-1.4° in latitude (latitudinal resolution is finer near the equator). There are 43 vertical layers. It is coupled to a dynamic-thermodynamic sea ice model. These models are used extensively for modern climate (Obase et al. 2017), paleoclimate (Abe-Ouchi et al. 2013Sherriff-Tadano et al. 2018), and future climate studies (Yamamoto et al. 2015). 50

MIROC4m-LPJ
We newly developed a vegetation coupled AOGCM MIROC4m-LPJ by introducing the Lund-Potsdam-Jenna Dynamical Global Vegetation Model (LPJ-DGVM; Sitch et al. 2003) into MIROC4m. The coupling method is the same as that used for MIROC-LPJ in previous studies (O'ishi and Abe-Ouchi 2009) which adopted a slab ocean instead of full ocean 55 model COCO. In the present study, we apply a direct coupling to avoid bias correction coupling between atmosphere and vegetation components. LPJ-DGVM predicts potential vegetation distribution which is translated to the classification used for MATSIRO annually by using monthly mean temperature, precipitation and cloud cover predicted by the atmosphere component of the GCM. A detailed description of the method can be found in O'ishi and Abe-Ouchi (2009). Another important modification is the introduction of a wetland scheme developed by Nitta et al. (2017). This scheme improves seasonality of 60 the hydrological behaviour over land. When snow melt occurs, a part of melt water does not directly runoff into river but stored https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License.
in an isolated storage. This stored water decays with timescale depending on the standard deviation of topography and the decayed amount of water is taken into account of land surface water and energy balance. The introduction of this scheme reduces summer warm bias over land in middle to high latitude. The model resolutions are the same as those of MIROC4m.

Settings
In the present study, three models are run with the same forcings and with the boundary conditions of the pre-industrial 80 (piControl in PMIP4/CMIP6, hereafter PI) and of the Last Interglacial (lig127k in PMIP4, hereafter 127k), as shown in Table   1 Vegetation distribution in MIROC4m is fixed to a present-day vegetation according to Ramankutty and Foley (1999); see  (Matthews 1983, Matthews 1984, 90 Hall et al. 2006 as the same as that in MATSIRO of MIROC5 (Watanabe et al. 2010) (see Figure 2) and as the same as that in the DECK piControl experiment (Hajima et al. 2019) using VISIT-e. However, as described above, VISIT-e predicts the leaf area index (LAI) which is referred by MATSIRO. In all experiments, the last 100 years, during which the climate has reached equilibrium, are used for analysis. https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License.

Temperature
Seasonally and annually averaged surface air temperature difference between the 127k from PI are shown in Figure 4. All models show the largest warming (>6K) in June-July-August (JJA) and global cooling in December-January-February (DJF) which corresponds to increased boreal summer solar irradiance and decreased boreal winter solar irradiance in the LIG, respectively. In September-October-November (SON), there is a switch to cooling over northern hemisphere land at middle 100 and high latitudes. On the other hand, warming (>3K) occurs across Antarctica in all models, reflecting the increase of solar irradiance at SH high latitudes in SON. In March-April-May (MAM), the global cooling is less than that in DJF and warming is seen at northern high latitudes.
The change in annually averaged surface air temperature is smaller than that of the seasonally averaged value because summer May-June-July-August (MJJA) averaged temperature difference is shown in Figure5, focusing on Greenland. MIROC4m and MIROC4m-LPJ show similar pattern of warming which is larger as the altitude is higher. MIROC-ES2L shows more homogeneous warming pattern over Greenland and intensity of warming is weaker than former two models.
Simulated temperature changes are compared with reconstructed values from proxies. In Figure 6, annual surface temperature 115 change is compared with land proxies by Turney and Jones (2010), hereafter referred to as TJ2010. MIROC4m-LPJ exhibits the largest warming among three models at northern high latitudes. However, the agreement with TJ2010 is not quantitative but just qualitative; annual warming is still underestimated. The other two models show cooling in northern high latitude land.
Greenland warming appears in all models, but is smaller than that of TJ2010. In Antarctica, MIROC-ES2L shows the same sign of temperature change as in TJ2010, but the intensity (at most +1K) is weaker than that of TJ2010. MIROC4m and 120 MIROC4m-LPJ show cooling rather than warming in Antarctica. In Figure 7, annual surface temperature change is compared with a newer reconstruction by Capron et al. (2017), hereafter referred to as C2017. All models show warming in Greenland, but only MIROC-ES2L reproduces Antarctica warming. The models underestimate warming at all sites in C2017 as the intensity is not reproduced. Figure 8 compares the simulated annual sea surface temperature (SST) change and TJ2010 ocean proxies. All three models 125 predict warm SST at northern high latitudes and cooling at low latitudes. This large-scale latitudinal pattern agrees with TJ2010, but some individual sites disagree in terms of sign. For example, the NH SST warming is more substantial in MIROC4m-LPJ https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License. and in MIROC-ES2L than in MIROC4m. The largest warming in the SH is seen in MIROC-ES2L. However, the intensity of SST changes in all models is far smaller than that of TJ2010. We also compare model annual SST difference with newer reconstructions by C2017 and Hoffman et al (2017), hereafter referred to as H2017, in Figure 9. MIROC4m underestimates 130 warming in the northern Atlantic Ocean and shows changes of opposite sign in the southern part of the Pacific and Indian Oceans. MIROC4m-LPJ shows larger warming in NH than in MIROC4m. However, as with MIROC4m, warming in SH is not simulated. MIROC-ES2L predicts better warming in the northern Atlantic Ocean than those of MIROC4m/MIROC4m-LPJ. MIROC-ES2L also predicts improved warming in SH which are partially consistent with proxies. Summer temperature change in the models are compared with that of reconstructions by C2017 and H2017 in Figure 10. Across the wide expanse 135 of the northern Atlantic Ocean, all models predict warming whose sign is consistent with that of reconstructions, except for some cooling sites. MIROC4m-LPJ predicts the largest warming (at most > +4K) while the other two models show a smaller intensity of warming (at most +3K) in the northern Atlantic Ocean. On the other hand, in SH, MIROC4m and MIROC4m-LPJ show cooling in contrast to the warming indicated by proxies. MIROC-ES2L show warming across much of the Southern Ocean, in contrast to the other two models. However, some sites still indicate an opposite sign to that of proxies. 140

Precipitation
The precipitation change between 127k and PI is indicated in Figure 11. In general, similar annually averaged precipitation changes are seen in all three models. Precipitation increase occurs in subtropical arid areas in the northern hemisphere, 145 especially in Africa. MIROC4m and MIROC4m-LPJ show almost the same precipitation change since they share the same atmosphere component. MIROC-ES2L shows a pattern slightly different to that of the MIROC4m models. Figure 12 shows JJA zonally averaged precipitation over Sahara (30W-20E). In all models, 127k summer precipitation shifts northward compared to PI. MIROC4m-LPJ shows the largest northward shift of precipitation. MIROC4m predicts slight northward expansion of vegetation compared with that of PI reflecting the precipitation increase. However, vegetation 150 distribution in Sahara does not show a significant northward expansion known as the "green Sahara".

Sea ice
In figure  models. The sea ice thickness is smaller in MIROC4m-LPJ than in MIROC4m although they show a similar distribution pattern. This is due to different equilibrium temperatures in MIROC4m and MIROC4m-LPJ. The PI equilibrium temperature in MIROC4m-LPJ is slightly higher (not shown) than that of MIROC4m due to the inclusion of a vegetation component. Figure   14 shows September NH sea ice thickness. As well as in March, MIROC-ES2L shows a different distribution of sea ice and 165 response in 127k compared to MIROC4m-based models. MIROC4m-LPJ shows the smallest sea ice thickness among the three models in both PI and 127k which corresponds to a warm Arctic Ocean in both PI and 127k. MIROC4m-LPJ predicts obviously less sea ice in September in PI compared to observation (e.g. HadISST; Rayner et al. 2003). March SH sea ice thickness is shown in Figure 15. In the Southern Ocean, all models show sea ice thickness to be thinner than that in the Arctic Ocean.
March sea ice shows the same characteristics in the two MIROC4m-based models. MIROC-ES2L shows a different pattern, 170 and the smallest amount of sea ice, compared with MIROC4m-based models. In all three models, March sea ice increases in 127k but the intensity differs depending on the model as well as NH. September SH sea ice thickness is shown in Figure 16.
As in March, there is a discrepancy between the sea ice thickness distributions of the MIROC4m-based models and MIROC-ES2L. Sea ice in the PI is smaller in MIROC-ES2L and the response of sea ice in 127k is also the smallest in that model.

Vegetation
We compare the vegetation distribution in all three models (Figure 3) regardless of different treatment (prescribed or predicted).

MIROC4m and MIROCES2L adopt a fixed vegetation distribution based on satellite data. Vegetation distribution in
MIROC4m is based on the classification of MATSIRO, translated from actual vegetation by Ramankutty and Foley (1999). 180 MIROC-ES2L is also fixed to satellite-based vegetation distribution which is translated from satellite data (Matthews 1983, Matthews 1984, Hall et al. 2006. These two vegetation maps show similar patterns of forest and grassland, although differ in the interpretation of classification such as C3/C4 or the boundary between forest and tundra. Only MIROC4m-LPJ predicts vegetation distribution in the present study. The 100-year averaged vegetation distribution in the PI shows characteristics common with the other two satellite-based distributions, except for the overestimation of forests (in boreal forest band and 185 African tropical forest) and underestimation of grassland (in African Savanna and central Eurasia). In the 127k simulation, vegetation drastically changes at northern high latitudes. Tundra is broadly replaced by boreal deciduous forest and almost disappears, reflecting the summer warming in the northern high-latitude land especially in Eastern Siberia and North America.
Forestation of tundra regions causes amplification of warming at Eastern Siberia and North America, especially in the snow melting season. This northward shift of boreal forest finally affects by +3K annually averaged temperature in Eastern Siberia 190 and North America compared with LIG warming without vegetation change. Grassland appears at the wide area of borealtemperate boundary in both Eurasia and North America due to less precipitation to support forests in 127k. This increase of grassland causes cooling at middle latitudes, especially in Eurasia. In Sahara, slight northward expansion of grassland is seen, but MIROC4m-LPJ did not reproduce so called "green Sahara". https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License.

Conclusions and Discussions 195
In the present study, we examined the LIG and PI simulations in accordance to the PMIP4 protocol by using three different versions of the MIROC AOGCM. These three models show basically the same response of temperature to the LIG boundary conditions, ie warming in boreal summer and cooling in boreal winter (Figure 4). However, the annually averaged temperature is different among the models. Only MIROC4m-LPJ predicts annual warming at NH high latitudes qualitatively consistent with proxy data such as TJ2010, C2017 and H2017, while the other two models show a cooling at NH high latitudes. 200 The vegetation change seen in MIROC4m-LPJ simulations is a reasonable response to temperature change induced by orbital parameter. The largest change is the northward shift of boreal forest and expansion of grassland in middle latitude. By comparing MIROC4m-LPJ and MIROC4m, we suggest that the vegetation feedback mechanism is necessary to explain the temperature change reconstructed by proxies since MIROC4m-LPJ predicts closer warming to reconstructions. By considering the overestimation of forest in PI, vegetation feedback may be still underestimated in MIROC4m-LPJ. The introduction of 205 dynamical vegetation in MIROC4m-LPJ appears to amplify the warming not only over land but also in the ocean at NH high latitudes. On the other hand, MIROC-ES2L, which partially introduces a vegetation effect through LAI prediction, does not show enough warming in LIG, but even shows annual cooling over land in the northern high latitude. This result indicates that vegetation feedback including change of vegetation distribution is necessary to predict warm climate both in past and in future.
Moreover, it is also pointed out that even MIROC4m-LPJ underestimate warming in spite of the large climate sensitivity (about 210 4 deg C) of AGCM of MIROC4m among GCMs (Meehl et al. 2007).
Compared to observation, MIROC-ES2L shows most realistic PI distribution of sea ice in the Arctic Ocean, owing to a new ice physics model. MIROC4m-LPJ predicts the smallest amount of sea ice in PI among three models in both March and September, because temperature in MIROC4m-LPJ is generally higher than that of MIROC4m over land due to the inclusion of dynamical vegetation. This warm bias reduces sea ice in the Arctic Ocean in the pre-industrial and thus inevitably affects 215 on the warming in LIG. To investigate the mechanisms in detail, we are planning further feedback analysis focusing on surface energy balance.

Code and Data availability
The code of MIROC-ES2L, MIROC4m, MIROC4m-LPJ are not publicly archived because of the copyright policy of MIROC 220 community. Readers are requested to contact the corresponding author if they wish to validate the model configurations of MIROC family models and conduct replication experiments. The source codes, required input data, and simulation results will be provided by the modeling community to which the author belongs. The output of the piControl and lig127k from MIROC-ES2L will be distributed and made freely available through the Earth System Grid Federation (ESGF). Details on the ESGF can be found on the website of the CMIP Panel (https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6). All experiments 225 done by MIROC4m and MIROC4m-LPJ will be available from following FTP server ( Figure 5: May-June-July-August averaged temperature difference (K) between 127k and PI in three models. The present-day calendar is applied to both PI and 127k. 445 Figure 6: Annual surface air temperature change between 127k and PI is compared with reconstruction by Turney and Jones (2010) https://doi.org/10.5194/cp-2019-172 Preprint. Discussion started: 17 January 2020 c Author(s) 2020. CC BY 4.0 License. Figure 7: As same as Figure 6 but compared with reconstruction by Capron et al. (2017).