Set-up of the PMIP 3 paleoclimate experiments conducted using an Earth system model , MIROC-ESM

Paleoclimate experiments using contemporary climate models are an effective measure to evaluate climate models. In recent years, Earth system models (ESMs) were developed to investigate carbon cycle climate feedbacks, as well as to project the future climate. Paleoclimate events can be suitable benchmarks to evaluate ESMs. The variation in aerosols associated with the volcanic eruptions provide a clear signal in forcing, which can be a good test to check the response of a climate model to the radiation changes. The variations in atmospheric CO 2 level or changes in ice sheet extent can be used for evaluation as well. Here we present implementations of the paleoclimate experiments proposed by the Coupled Model Intercomparison Project phase 5/Paleoclimate Modelling Intercomparison Project phase 3 (CMIP5/PMIP3) using MIROC-ESM, an ESM based on the global climate model MIROC (Model for Interdisciplinary Research on Climate). In this paper, experimental settings and spin-up procedures of the mid-Holocene, the Last Glacial Maximum, and the Last Millennium experiments are explained. The first two experiments are time slice experiments and the last one is a transient experiment. The complexity of the model requires various steps to correctly configure the experiments. Several basic outputs are also shown.


Introduction
For better future climate projection, evaluating climate models using experiments over different time scales is gaining recognition.As an experimental design, paleoclimate simulations provide unique opportunities to test models.Recent periods that are climatologically quite different to present-day conditions are targeted for the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al., 2011)/Paleoclimate Modelling Intercomparison Project phase 3 (PMIP3).
Here we describe three paleoclimate experiments using MIROC-ESM -an Earth system model (ESM) based on MIROC (Model for Interdisciplinary Research on Climate) -which are conducted in the CMIP5/PMIP3 framework and for the Intergovernmental Panel on Climate Change 5th Assessment Report (IPCC AR5).
In this paper, we consider ESMs as coupled climate models with biogeochemical components which are capable of treating the carbon cycle online together with other climate components: in other words, climate models which can predict the atmospheric CO 2 concentration.Among such ESMs, the models of higher complexity, consisting of atmosphereocean general circulation models (AOGCMs) coupled with biogeochemical components, are now available for paleoclimate simulations.These models allow evaluation of more components compared to the previous AGCM work in PMIP phase 1, or AOGCM inter-comparison in PMIP phase For IPCC AR5, three paleoclimate experiments were officially listed in the CMIP5/PMIP3 project (Taylor et al., 2009;Otto-Bliesner et al., 2009): mid-Holocene (6000 yr before present time slice, hereafter 6 ka), Last Glacial Maximum (21 000 yr before present time slice, hereafter LGM), and Last Millennium (850-1850 AD transient, hereafter LM).We briefly review the scientific questions that draw the attention of the paleoclimate community.
The pronounced characteristics in the climate reconstructed for 6 ka are the enhancement of African and Asian summer monsoons and the associated pronounced vegetation changes over the Sahara desert (Hoelzmann et al., 1998;Jolly et al., 1998;Peyron et al., 2006).The enhancement and northward shift of the summer African and Asian monsoons were reported in PMIP phase 1 (Joussaume et al., 1999) and PMIP phase 2 (Braconnot et al., 2007).These changes were consistent in sign with paleoclimate proxy records, mostly from pollen (Kohfeld and Harrison, 2000;Prentice et al., 2000;Yu et al., 2000;Harrison et al., 2001Harrison et al., , 2003;;Bigelow et al., 2003;Pickett et al., 2004;Bartlein et al., 2010), but most of the models failed to simulate sufficiently high Holocene precipitation to maintain a vegetated Sahara (e.g.Doherty et al., 2000;Chikira et al., 2006).Much modelling work has investigated the vegetation feedback effect on the African monsoon change in 6 ka, and most have described the presence of positive feedback between vegetation and precipitation (Ganopolski et al., 1998;Doherty et al., 2000;Levis et al., 2004).Several studies have investigated the role of the ocean on the summer monsoons using AOGCMs, and most suggested that there is further enhancement of the African monsoon when the ocean is coupled.As for the Asian monsoon, it is still enhanced compared to 0 ka, but the signal is weaker due to the ocean feedback (Liu et al., 2003;Zhao et al., 2005Zhao et al., , 2007;;Ohgaito andAbe-Ouchi, 2007, 2009;Zhao and Harrison, 2011).Dallmeyer et al. (2010) used an ESM to investigate the effect of the ocean and vegetation component on the total change of the Asian monsoon, and found that the ocean plays a major role in the enhancement of precipitation.This demonstrates that experiments with ESMs coupled with dynamic global vegetation models provide an opportunity to evaluate the comprehensive influence on Asian and African monsoon change resulting from dynamic vegetation and oceans.
The LGM is characterized by a very cold and dry climate.The Laurentide and Fennoscandian ice sheets cover North America and Northern Europe (e.g.Denton and Hughes, 1981), atmospheric greenhouse gas (GHG) levels are reduced, and dust transportation is enhanced (e.g.Barnola et al., 1987;Petit et al., 1981).Land vegetation responds as well: subtropical deserts expand and forests generally regress (Colinvaux et al., 1996(Colinvaux et al., , 2000;;Marchant et al., 2009).At high latitudes, boreal forests migrate southward, and are replaced by tundra and grassland (Prentice et al., 2000;Tarasov et al., 2000;Ray and Adams, 2001;Harrison and Prentice, 2003).Due to these significant differences from today, various aspects of the climate can be investigated.For example, attempts to constrain climate sensitivity using LGM cooling (Annan et al., 2005;Schmittner et al., 2011;Hargreaves et al., 2012), and the mechanism of the recorded weakening of the Atlantic meridional overturning circulation (AMOC) and the associated cooling in the North Atlantic Ocean (Weber et al., 2007;Murakami et al., 2008), have been major subjects since PMIP2.As for the atmospheric circulation, the southward shift of the Intertropical Convergence Zone (ITCZ) and changes in the Southern Hemisphere westerlies (Rojas et al., 2009) have also been studied.Temperature changes have been discussed in comparison with archives of proxy records (Otto-Bliesner et al., 2009;MARGO project members, 2009;Bartlein et al., 2010).With respect to the carbon cycle, the glacial-interglacial atmospheric CO 2 variations remain an open issue (Kohfeld and Ridgwell, 2009).LGM climate experiments using ESM with an interactive carbon cycle will contribute to the understanding of the mechanisms controlling CO 2 variations.
The LM includes two key climatic epochs: the "Medieval Climate Anomaly" (MCA, ca.1000-1200 AD) and the "Little Ice Age" (LIA, ca.1500-1850 AD).Recent progress on proxy-based climate reconstruction has slowly revealed the spatial and temporal extent of those climatic epochs (e.g.Trouet et al., 2009;Mann et al., 2009), which allows better comparison between data and GCM simulations.Such comparison is useful in verifying the ability of a climate model, and assists the understanding of the mechanisms contributing to the climate variability.More specifically, the need to distinguish anthropogenic climate change from natural climate variability at centennial to millennial time scales is one of the key issues of climatology, addressing the social demand for climate research.Several very large volcanic eruptions occurred during this period, causing global cooling as a response to the radiation changes.Combined with the reconstructions, they serve as ideal test cases to investigate the response of climate models.Volcanic activity may also play a key role in reproducing LIA cooling, by decadally paced volcanic activity starting around the 13th century (Zhong et al., 2011).LM experiments with fully coupled ESMs allow discussion of the evolution of the carbon cycle interacting with the climate variation over the last 1200 yr.
LM simulations appearing in the IPCC AR4 (IPCC Fourth Assessment Report, 2007) were performed mostly using Earth system models of intermediate complexity (EMICs) (e.g.Gerber et al., 2003;Goosse, 2005).Following some work with coupled models without a carbon cycle (e.g.Gonzáles-Rouco et al., 2003;Ammann et al., 2007;Servonnat et al., 2010), a LM experiment using an ESM has been reported recently (Jungclaus et al., 2010).Our study uses a similar approach in evaluating the strength of natural variability using up-to-date forcing reconstructions.This paper presents the technical aspects of CMIP5/PMIP3 paleoclimate experiments using MIROC-ESM, including model settings, choice of boundary conditions, special treatment for model spin-up, and the initial data preparation.It is not our intention to present scientific results in detail here, as these will be discussed in subsequent papers.In our preliminary results, we show the seasonal changes for 6 ka, whereas annual means for 21 ka are presented.The forcing difference for 6 ka is principally the orbital configuration, which causes the change in seasonal pattern of the incoming shortwave radiation.As for 21 ka, the different equilibrium state is achieved responding to the continental-scale topography and global albedo change, which should be represented by the annual means.
The model and its general characteristics are described in Sect.2, specific settings for each paleoclimate experiment are described in Sect. 3 together with preliminary results, and concluding remarks follow in Sect. 4.

General characteristics of the model
MIROC-ESM is an ESM developed at the Japan Agency for Marine-Earth Science and Technology in collaboration with the Center for Climate System Research (CCSR), University of Tokyo, and the National Institute for Environmental Studies based on the MIROC AOGCM (K-1 model developers, 2004;Nozawa et al., 2007).Based on MIROC3.2, which is the previous version of MIROC used for the CMIP3 experiments, two versions are now available as AOGCM: the bug-fixed version MIROC4 and the new generation model MIROC5.MIROC-ESM is another branch of the MIROC series based on MIROC4, in which ecosystem components are introduced as an ESM.MIROC-ESM consists of an atmospheric general circulation model (MIROC-AGCM 2010), including an online aerosol component (SPRINTARS 5.00), an ocean GCM with sea ice component (COCO 3.4), a land surface model (MATSIRO), and terrestrial and marine ecosystem components.These are interactively coupled in MIROC-ESM, as illustrated in Fig. 1.These atmosphere, ocean, and land surface components, as well as a river routine, are coupled by a flux coupler (K-1 Model Developers, 2004).As an ESM, MIROC-ESM has a newly introduced atmospheric chemistry component (CHASER 4.1) and carbon cycle components for the land and ocean ecosystems.The carbon cycle is calculated by a marine biogeochemical component coupled with COCO (CCSR Ocean COmponent model) and a terrestrial ecosystem component dealing with dynamic vegetation (SEIB-DGVM) coupled with MAT-SIRO.Watanabe et al. (2011) gives detailed descriptions of each component and the model's performance against observed 20th-century climate variations.
The atmosphere and land components have a spatial resolution of T42 (equivalent grid interval is approximately 2.8 • in latitude and longitude).The atmosphere has 80 vertical layers between the surface and the stratosphere up to about 0.003 hPa (Watanabe et al., 2008).Due to the limited computational resources, part of the CMIP5 experiments were performed with the full version of MIROC-ESM, to which atmospheric chemistry component (CHASER) is coupled.Other experiments, including the paleoclimate experiments described in this paper, were performed without CHASER.
The ocean component COCO has a resolution of 1.4 • in longitude and 0.56-1.4• in latitude in the horizontal and 44 levels in the vertical (K-1 Model Developers, 2004;Hasumi, 2000).No correction is applied in exchanging heat, water, and momentum flux between the atmosphere and the ocean.The sea ice model is based on a two-category thickness representation, has zero-layer thermodynamics (Semtner, 1976), and dynamics with elastic-viscous-plastic rheology (Hunke and Dukowicz, 1997).
The land surface component MATSIRO (Minimal Advanced Treatments of Surface Interaction and RunOff) consists of a single layer canopy, three layers of snow, and six layers of soil to a depth of 14 m (Takata et al., 2003).It has 11 types of land cover and soil (represented as land index; see Fig. 2), which classifies the ground surface conditions and provides the table of parameters for computing physical land processes.It is coupled to a river-routing model, TRIP (Total Runoff Integrating Pathways) (Oki and Sud, 1998), for calculating river discharge.The ageing effect on snow albedo  and the effects of dirt in snow have been taken into account (Yang et al., 1997;Aoki et al., 2006).The dirt concentration in snow is calculated from the deposition fluxes of dust and soot calculated in the aerosol module SPRINTARS (Spectral Radiation-Transport Model for Aerosol Species) (Takemura et al., 2000(Takemura et al., , 2002;;Takemura, 2005).
MATSIRO is coupled with the terrestrial ecosystem model SEIB-DGVM (Spatially Explicit Individual-Based Dynamic Global Vegetation Model) (Sato et al., 2007;Ise et al., 2009).SEIB-DGVM adopts a scheme that explicitly captures light competition among trees, avoiding parameterized competition.Land vegetation is classified into 13 plant functional types (PFTs), and each PFT has different ecophysiological parameters and allometry relationships.This results in differential growth patterns and competition among PFTs under the environmental conditions in each grid cell.Leaf area index (LAI) is predicted at daily time steps by computing plant physiological processes, such as photosynthesis, respiration, growth, and mortality, based on the climatic conditions simulated by MIROC-ESM.The predicted LAI is then used for calculation of the biogeophysical processes in MAT-SIRO: the processes for radiation transfer including surface albedo, interception of precipitation, and energy transfer by latent/sensible heat.In the current configuration, MATSIRO considers the dynamic variations in LAI, while the land index in MATSIRO is fixed, and neglects the PFT predicted by SEIB-DGVM.Thus, the terrestrial ecosystem in MIROC-ESM affects climate only through the carbon cycle and the LAI feedback.
The oceanic biogeochemical model is based on a simplified version of a nutrient-phytoplankton-zooplanktondetritus (NPZD) ecosystem model (Oschlies, 2001;Oschlies and Garc ¸on, 1999).The oceanic biogeochemical tracers include nitrate, phytoplankton, zooplankton, detritus, calcium carbonate, calcium, alkalinity, and dissolved inorganic carbon (DIC).A constant Redfield ratio (C/N = 6.625) is used to estimate the carbon and calcium flow.The sea-air CO 2 flux is calculated by multiplying the difference of oceanatmosphere CO 2 partial pressures by the ocean gas solubility.These model components are all activated in our paleoclimate experiments which will be described in the next section.

Experimental setup
The basic settings of the experiments are common with other simulations in CMIP5 by MIROC-ESM (Watanabe et al., 2011).Here, the pre-industrial control experiment (PI) is firstly described as a counterpart of the following paleoclimate simulations.PI is a time slice experiment forced by the climatic conditions at 1850 AD.The greenhouse gas levels (GHGs) are set to 285 ppm for CO 2 , 0.3 ppm for N 2 O and 1.7 ppm for CH 4 (Table 1).The land index for MATSIRO is shown in Fig. 2a.A plausible distribution in 1850 is prepared based on modern vegetation types and historical land use data.The orbital parameters are set to 0.01672 for the eccentricity, 23.45 • for the obliquity, and 102.04 • for the angular precession.Although these orbital parameters have a considerable effect on longer time scales, they have had little impact on annual mean insolation on time scales shorter than the millennium.According to IPCC AR4, summer insolation decreased by 0.33 W m −2 at 45 • N over the millennium, winter insolation increased by 0.83 W m −2 (Goosse, 2005), and the magnitude of the mean seasonal cycle of insolation in the Northern Hemisphere decreased by 0.4 W m −2 .

Initialization of the PI simulation
Initial data for the PI experiment was prepared in stages.The procedure is presented in detail by Watanabe et al. (2011), the description paper of MIROC-ESM.Here we outline the procedure briefly in Fig. 3a.Starting with the "lowtop" version of MIROC-ESM with 20 vertical layers in the were separately performed beforehand and combined.A 50 yr integration followed with an atmosphere with 80 vertical layers.In the course of the spin-up runs, representative states and fluxes in the physical climate (surface air temperatures, radiation fluxes at the top of the atmosphere, strength of the thermohaline circulation, sea ice extent) and carbon cycle components (soil and vegetation carbon storage) were monitored.The spin-up run was continued until linear trends in the last 50 yr of these quantities became insignificant.Although the climatic trends were thought to be insignificant, a slight warming drift (0.1 • C century −1 in global average surface air temperature) still existed in the PI experiment, which was overlooked during the spin-up procedure.It is considered to be due to the slow reduction in sea ice extent (of 3 % century −1 ) in the Southern Ocean, which caused a similar warming drift in the LM experiment (Sect.3.4).Because of this drift, additional care is needed in analysis.After these spin-up procedures, a 630 yr integration was carried out as the PI experiment.The last 100 yr of the PI run are used as a reference for the analyses of the time slice experiments (6 ka and LGM).The trend is +0.05 • C in this period.

Model performance
In this section, the climate of the PI experiment is briefly presented comparing with the modern climate.Although the PI simulation is not necessarily comparable with modern-day observations, the aim of this section is to describe the characteristics of the PI climate, which is commonly used as a baseline to compare the paleoclimate experiments with proxy data.A detailed comparison of a present-day simulation with modern observations was already presented by Watanabe et al. (2011)  The sea surface temperature (SST) is shown in Fig. 4 in comparison with the World Ocean Atlas temperature data (WOA, 1998;Levitus et al., 1998).PI simulated reasonable SST distribution, though biases appear in the Pacific Convergence Zone, the Kuroshio region, and the eastern Equatorial Pacific.In the global average, slightly cooler SST in PI is reasonable considering the 20th-century warming.The distributions of simulated precipitation in comparison with observational data (Xie and Arkin, 1996) for boreal summer (JJAS) and winter (DJF) are shown in Fig. 5.The precipitation is underestimated along the South Pacific Convergence Zone (SPCZ) and over Central America, whereas it is overestimated over the Maritime Continent and the northwestern Indian Ocean.These shortcomings are similar to those in our previous model (MIROC3.2) because these two models have almost the same atmospheric physics components.For comparison with paleoclimate, reasonable representation of the ITCZ and precipitation distribution of the monsoon area is encouraging.
Further detailed comparison of observational data with the 20th-century part of the LM experiment would be beyond the scope of this paper, which aims at the description of the experimental settings.

Forcing data and boundary conditions
The 6 ka time slice experiment was performed following the PMIP3 protocol (Table 1).The boundary conditions are kept unchanged from the PI run except for the orbital parameters (eccentricity: 0.018682; angular precession: 0.87 • ; obliquity: 24.105 • ) and GHG concentrations (CO 2 : 280.0 ppm; N 2 O: 0.65 ppm; CH 4 : 0.27 ppm).No volcanic activity was considered.
Since 6 ka GHG levels are not drastically different from PI condition (except for methane), this experiment essentially tests the climate response to a change in seasonal incoming shortwave radiation due to different orbital parameters (Fig. 6a).Incoming shortwave radiation is increased in boreal summer and decreased in austral summer.As for annual radiation, there is an increase in the high latitude, and a decrease in the tropics.Through these changes, seasonality is enhanced in the Northern Hemisphere and weakened in the Southern Hemisphere, especially at mid-latitudes.

Preparation of the initial data
Figure 3b shows the procedure to set up the 6 ka experiment.It has been branched off year 250 of the PI run by replacing GHG levels and the orbital parameters with 6 ka values.
As the land surface condition in the PI experiment is based on the condition in year 1850, some grid cells are classified as cropland.To perform the 6 ka experiment, we need to replace such grid cells with the natural vegetation and to spin up the terrestrial ecosystem to a steady state without land use for cultivation.To shorten the spin-up time, the initialization  was conducted with the offline terrestrial ecosystem model that is used in MIROC-ESM.Firstly, 100 yr of integration was performed with the 6 ka GHG concentrations and orbital forcing by MIROC-ESM.Using the last 25 yr of this output, we ran the offline model for 2000 simulation years with the cycling 25 yr climate forcing.This procedure is a similar approach to the one described in the protocol of the Coupled Climate Carbon Cycle Model Intercomparison Project (C4MIP, Fung et al., 2000;Cox et al., 2002).After this procedure, the quantities of terrestrial ecosystems were merged back into the restart data of MIROC-ESM, followed by a spin-up of 180 simulation years with all components coupled (Fig. 3b).

Results
The results of the 6 ka experiment are also influenced by the warming drift identified in the PI experiment.In order to minimize the effect of drift, the following analyses are based on the anomalies between 6 ka and PI, and the integration time from the point of initialization of 6 ka (i.e."branch" time) was kept the same for both experiments.
The changes in simulated temperature 6 ka PI for JJAS and DJF are shown in Fig. 7. Pronounced changes are the warming over most of the boreal continents by 1-3 • C in JJAS in response to the strengthened solar radiation in the boreal summer.The tropical oceans are slightly cooled despite the increased incoming radiation, which is mainly due to the 2-3-month lags in ocean response (Braconnot et al., 2000).In DJF the cooling over the continents and the ocean are simulated as a response to the change in incoming radiation.The continents are cooled stronger by 1-4 • C, and the ocean is cooled by 0.5-1 • C. The result is consistent with the PMIP2 multi-model ensemble (Braconnot et al., 2007).
The changes in simulated precipitation are shown in Fig. 8.In JJAS, increases of precipitation over the Sahel and northern India by 1-3 mm day −1 and the reduction of precipitation south of these areas in similar magnitude are simulated.These precipitation changes and the circulation changes suggest an enhancement of the African and Asian monsoon.However, while these enhanced monsoons are consistent with the proxy records, the amount of precipitation over most of the Sahara Desert is not enough to explain the proxy records.In monitoring the daily predicted PFTs over the Sahara, SEIB-DGVM predicts non-desert vegetation more often than in the PI case, though non-desert PFTs are still much less frequent than desert vegetation.The monsoon activities for the boreal summer are investigated in another work (Ohgaito et al., 2012).In the Southern Hemisphere, a weakening of precipitation over the continents by 0-2 mm day −1 is simulated in DJF, which suggests the weakening of monsoon activities.

Forcing data and boundary conditions
The LGM experiment is an equilibrium experiment using the boundary conditions defined by the PMIP3 protocol (Table 1) except for the salinity setting.The PMIP3 protocol recommends increasing the mean ocean salinity by 1 PSU everywhere at the beginning of the simulation.In our experiment, however, the salinity setting was unchanged from the PI condition.A salinity-modified experiment remains to be done in the near future.The change in seawater density leads to change in the mixed layer, which may cause a change in airsea interaction.Atmospheric CO 2 /N 2 O/CH 4 concentrations, as well as orbital parameters, are set to the values shown in Table 1.The shortwave solar radiation deviation at the top of the atmosphere shown in Fig. 6b is much smaller compared to the 6 ka orbital change.
Another major difference in the LGM boundary condition is the topography.The Laurentide and Fennoscandian ice sheets covered wide areas of North America and the northern half of Europe.In the ESM, ice sheets are expressed as mountain ranges with albedo of ice (i.e.flow is not considered).Three ice sheet reconstructions are currently available: ICE-6G (Peltier, 2009;Argus and Peltier, 2009;Peltier and Drummond, 2008;Peltier et al., 2013), ANU Ice Model (Lambeck and Chappell, 2001;Lambeck et al., 2002Lambeck et al., , 2003)), and GLAC-1 (Tarasov andPeltier, 2002, 2003).Since those reconstructions vary in topography, the settings of our experiment follow the "PMIP standard ice sheet" recommended by the PMIP3 committee, which averages the three reconstructions (PMIP3, 2010; Abe-Ouchi et al., 2013).The data provided by PMIP3 were regridded to T42 resolution and transformed to the spectral wave number space in order to adjust to the atmosphere and land component.The orographic difference for LGM-PI is shown in Fig. 9.  Sea level was lowered by about 90 m, which closes the Bering Strait, and joins together the Indonesian maritime continent.Sea level change is set to 90 mm to maintain the consistency with the PMIP3 topography reconstruction (i.e.ice sheet topography).The total mass of the atmosphere was not adjusted.The vegetation boundary condition of the LGM simulation for MATSIRO was created based on the PI vegetation distribution.At first, the vegetation types in the ice sheet area given by the PMIP3 protocol are changed to ice sheet, and the rest of the land areas are basically not changed from PI.However, due to the lower sea level in the glacial climate, the area of the continental shelves should have some vegetation type specified, and the same vegetation types from neighbouring grid cells in zonal direction are given (Fig. 2b).

Preparation of the initial data
The initial data was prepared following the procedure shown in Fig. 3c.The initial condition of the ocean physical field is inherited from the result of the LGM experiment following the PMIP2 protocol with an earlier version of MIROC, MIROC4m (a bug-fixed version of MIROC3.2 (medres)), having the same resolution.This experiment has been integrated for 1900 yr under PMIP2 LGM conditions and is regarded as being in equilibrium.Since the PMIP3 LGM ice sheet topography differs from the PMIP2 version, a spin-up for the ocean is necessary.
As the atmospheric initial condition, the PI experiment by MIROC-ESM is used, as is the case for the other two experiments.Since the atmosphere of MIROC-ESM has a different vertical coordinate system from previous versions, it is difficult to initialize the model by applying the initial condition of these old versions.Thus the atmospheric part taken from the 250th year of the PI run is combined with the above-mentioned ocean part.Then we conduct the spin-up for 150 model years using the coupled model, with the orbital parameters, GHG levels, topography, and the land-sea mask of LGM conditions.Since the adjustment time of the atmospheric physical field is considerably shorter than that of the  ocean, the atmospheric initial condition does not affect the length of the spin-up.
Meanwhile, special care is needed to start an integration under the new topography, which has the large ice sheets in the Northern Hemisphere.The surface pressure field must be adjusted to the change in surface elevation over the continents.As MIROC-ESM has high vertical resolution in the stratosphere, the model was sensitive to the topography change, generating gravity waves and breaking the Courant-Friedrichs-Lewy condition by generating unrealistic high wind velocities.As recommended in the PMIP3 protocol, we changed the surface elevation gradually (in two stages, integrated one model year for adjustment), and also adjusted the initial pressure field to the LGM surface elevation when starting the integration under the new topography.
Since the atmospheric CO 2 concentration in the LGM experiment needs to be reduced by ca. 100 ppm compared to the PI condition, the carbon budget requires long integration to reach a steady state.After the physical field of ESM (i.e.atmosphere-ocean-land) achieved a quasi-equilibrium state, land and ocean biogeochemical components were run sepa-   climate to obtain a steady state.During this process, the land use (cropland) is also reset to the potential natural vegetation.Distribution of the ice sheets was taken into account for suppressing the new vegetation and accumulation of soil carbon over the ice grid cells.
The spin-up of the oceanic ecosystem was performed in a similar manner.Starting from the PI condition, the offline model (Chikamoto et al., 2012) was integrated for 3500 model years with the LGM climatology (i.e.not cyclic forcing).The partial pressure of the atmospheric CO 2 was set to 185 ppm, to calculate air-sea CO 2 gas exchange.This setting forces the near-surface oceanic pCO 2 to be 189 ppm and achieved a new equilibrium of marine carbon distribution.It results in the removal of the oceanic CO 2 , which is equivalent to an increase of ca.550 ppm CO 2 in the atmosphere.
After these offline experiments had achieved equilibrium, variables were merged back into the experiment using MIROC-ESM, and the model was integrated over 200 yr to obtain an equilibrium state for the whole system under the LGM condition.The final 100 model years are used for analysis.

Results
The change in simulated mean annual temperature for LGM-PI is shown in Fig. 10.More than 20 • C cooling over the Laurentide and Fennoscandian ice sheets is seen.The rest of the globe cooled much less.Tropical cooling ranges between 1 and 3 • C, which is consistent with the proxy records (MARGO Project Members, 2009).Comparing the temperature response over land to the ocean, the continents cooled more than the ocean as in earlier work (e.g.Braconnot et al., 2007Braconnot et al., , 2012;;Laîné et al., 2009).This relationship quantitatively agrees with other climate model outputs and proxy reconstructions (Schmidt et al., 2012).The changes in simulated annual precipitation are presented in Fig. 11.The precipitation was reduced for most of the simulated area, suggesting that a drier climate is associated with the cooling.The exception is the increased precipitation along the SPCZ, which is consistent with earlier work (e.g.Toracinta et al., 2004;Yin and Battisti, 2001).Toracinta et al. (2004) suggested that the difference in surface temperature gradients plays a major role in enhancing the SPCZ during LGM.
Although some proxy records suggest the weakening of the AMOC in LGM (McManus et al., 2004), the LGM experiment simulated ca.32 Sv, which is about 15 Sv stronger than that of PI (AMOC was defined by the maximum value of the stream function of the Atlantic meridional circulation at 30 • N).
LGM AMOC shows the stronger variation in amplitude (±3-5 Sv), with longer periods (5-8 yr), while PI AMOC shows weaker variation (±1-2 Sv, 3-4 yr).About half of the modelling studies have had a problem representing the weaker AMOC as suggested by the proxies (Weber et al., 2007).Oka et al. (2012) suggested that slight differences in surface cooling or wind stress forcing could lead to the different response of the AMOC.Further investigations are required to understand its behaviour.
The response of the carbon cycle and predicted vegetation in SEIB-DGVM are investigated in another work.

LM
The LM experiment is a transient experiment using timedependent boundary data, unlike the other experiments described in this paper.In addition, the experiment was carried out as a "CO 2 prognostic" run, so that the atmospheric CO 2 concentration was calculated by the model itself.It requires the additional procedure in preparing the data for atmospheric CO 2 concentration.In this section, the preparation/selection of these data is described.

Forcing data and boundary conditions
The forcing data and boundary conditions for the LM experiment include (a) orbital forcing, (b) solar forcing, (c) volcanic forcing, (d) non-CO 2 GHGs, and (e) cultivated land area (land use).In Fig. 12, the time series of volcanic and solar forcing data are shown (panel a and b).
Orbital forcing is based on a provided table of parameters (Berger, 1978), and the value is updated every year on 1 January.It affects the seasonality of incoming shortwave radiation, while total intensity is given by solar forcing data.As for solar forcing, multiple reconstructions of annual total solar irradiance (TSI) are available, and five reconstructions with different coverage periods are recommended in the PMIP3 protocol (Schmidt et al., 2012).Amongst these, we used Wang et al. (2005) (hereafter, WLS) for 1610-2000 AD and Delaygue and Bard (2010) (hereafter, DB) for 850-1609 AD (Fig. 12b).WLS is based on a spectral reconstruc-tion and on a flux transport model using the observed sunspot record, while DB is based on an Antarctica stack of 10 Be records.To splice WLS to DB, the amplitude of DB is scaled to the modern-to-Maunder Minimum TSI in the WLS reconstructions.The amplitude of high-frequency components (i.e.ca.11 yr period variation) is dependent on the solar activity.This resulted in very small amplitudes of high-frequency components during the 15th and 17th centuries.The radiation scheme in MIROC-ESM considers 29 spectral bands, in which spectrally resolved data provided by WLS can be used directly.The TSI data provided by DB are distributed consistently so that the spectrally integrated TSI agrees with the reconstruction.No enhanced variations in the UV region proposed by Lean et al. (2000) and used by Shindell et al. (2001) were applied.The parameterization for solar-derived ozone variations suggested by the PMIP3 protocol was not used either, to maintain the consistency with the already-conducted 20th-century experiments.
Two reconstructions of the volcanic forcing, Gao et al. (2008) and Crowley et al. (2008) as well as Crowley and Unterman (2012), are provided by the PMIP3 protocol, as time series of aerosol optical depth (AOD) variation.Both datasets are based on polar ice cores, but differ in their selection of ice cores, and the Gao data show stronger forcing in general.AOD estimates are based on a correlation between sulphate in the Antarctic ice cores and satellite AOD data, which were calibrated with the 1991 eruptions of Mt Pinatubo and Cerro Hudson (Sato et al., 1993).Aerosol size estimation is based on Pinto et al. (1989).In our experiment, Crowley data were used (Fig. 12a).The volcanic forcing is calculated using the AOD at 0.55 µm and the effective radius calculated in MIROC-ESM for scaling.The data are provided at four latitude bands of equal area.After the year 1850, Sato et al. (1993, updated) is used for AOD data.
CO 2 , CH 4 and N 2 O are considered as GHG forcing in the PMIP3 protocol, and a data table is provided (Joos, 2007;Schmidt et al., 2011).These data are derived from high-resolution ice cores in Antarctica (Battle et al., 1996;Etheridge et al., 1996Etheridge et al., , 1998;;Ferretti et al., 2005;MacFarling Meure et al., 2006;Flückiger et al., 1999Flückiger et al., , 2002;;Machida et al., 1995), connected to the instrumental data (Hansen and Sato, 2004) using a spline fit, with some adjustment (Joos and Spahni, 2008).For the atmospheric CO 2 concentration, however, we decided to perform LM as a prognostic CO 2 experiment in which CO 2 is calculated online (i.e.CO 2 concentration predicted by the carbon cycle component is used for radiation transfer calculation), aiming to discuss the response of carbon cycle to the solar and volcanic forcings during LM.Since the majority of CMIP5/PMIP3 LM runs are "CO 2 given", this was an unusual approach among CMIP5/PMIP3 experiments.The model predicts CO 2 using emission data.Anthropogenic emissions are considered only after 1850 AD, while emissions are assumed to be zero for the 850-1850 AD period.CH 4 and N 2 O are given as per the PMIP3 data table.
As a difference from the PMIP3 protocol, the cultivated land area was assumed to be unchanged.The land use for the 1850 condition was applied throughout the experiment including the spin-up period.

Preparation of initial data
As is the case with the 6 ka experiment, pre-industrial physical fields (i.e. from PI experiment) of the atmosphere and ocean were used for the initial condition of LM spin-up.Since the forcing conditions (orbital parameters, solar irradiance and GHG concentrations) are quite similar between the pre-industrial period and years around 850 AD, PI can be regarded as an initial condition for the 850 AD simulation.
Figure 3d shows the procedure of the spin-up.The 250th year of the PI experiment was used as the initial state, then 56 yr of integration were completed under the year-850 condition, with fixed GHG (including CO 2 ) concentration.Then CO 2 concentration was set to be free, so that the modelcalculated value is used in the radiation process of the atmosphere component.The initial value of the CO 2 concentration for the "CO 2 free" experiment was reset from 298.71 ppm to the level of 279.3 ppm (i.e.value at year 850; Joos, 2007), as the transient simulation was started.This resetting procedure was performed to maintain consistency with the reconstructed value, and also to avoid unnecessary CO 2 feedback.We adopted this procedure from the protocol of C4MIP (Fung et al., 2000;Cox et al., 2002).It breaks the conservation of the total carbon in a precise sense, but the effect is small enough to be neglected in evaluating the total budget.

Results and time series
Figure 12 shows the time series of forcing data and the simulated results.The volcanic and the solar forcings are shown in Fig. 12a and b, respectively, and the simulated annual mean surface air temperature in the Northern Hemisphere is shown in Fig. 12c.The warming drift in temperature is visible during the period until the 19th century (ca. 9 × 10 −2 • C 100 yr −1 ), while the warming during the 20th century is more pronounced (ca.0.5 • C 100 yr −1 ).The drift is also visible in the PI experiment, which is caused by a decreasing trend in the Southern Ocean sea ice.It is likely due to an insufficiency of the initialization.Removing the trend is necessary to discuss the long-term climate variations.
Figure 12e shows the same temperature time series as Fig. 12c but with the linear trend of warming drift observed in the PI experiment subtracted (0.094 • C 100 yr −1 ).The air temperature represents anomalies ( • C) from the 1971-to-2000 mean, and is averaged over the Northern Hemisphere for the comparison with proxy records.The slope of the trend was calculated from linear regression of the PI experiment results for the period between the year 250 and year 630 (i.e. after the "branching").The blue shade in Fig. 12e represents the range of temperature reconstruction between the 20th and 80th percentile by Frank et al. (2010).
Considering that the proxy data cannot reflect the highfrequency variations, the time series of the simulated temperature roughly follows the reconstructions for most of the period (Fig. 12e).However, contrary to the good agreement after the year 1700, some discrepancies can be seen in the earlier periods.For example, the simulated values do not show a clear indication of the LIA, while the reconstruction shows relatively low temperature in the 16th to 17th centuries.The mismatch is also large during the 12th century, where the simulated temperature tends to be higher than the reconstruction.The model also seems to overestimate the responses to the very large volcanic events (e.g.year 1229,1258,1809,1816), but this could be due to the limitation of detection capability of the reconstruction method (Mann et al., 2012).
Comparing the simulated temperature anomaly with the forcing data (panel a and b in Fig. 12), it is clear that volcanic aerosols have a strong influence on the surface air temperature, while the solar irradiance has minor effects.Under the current set-up of the model, LM climate variations should depend mainly on those two types of forcing.The faint LIA is considered to be due to the relatively constant solar forcing.
The simulated atmospheric CO 2 concentration (Fig. 12d, black line, shown with the reconstructed value in blue) is stable in the model over the pre-industrial period, while it shows a sharp increase after 1850, the period during which anthropogenic emissions are considered.It shows about 400 ppm near the surface in the year 2001, while the observed value in Mauna Loa was 370.13 ppm (Tans and Keeling, 2011) in the same year.Considering that the simulated value at year 1850 was 289.5 ppm, which is ca. 5 ppm higher than the reconstruction, the model still shows a tendency of slight overestimation.This is partly attributable to the warming drift of the global temperature.

Conclusions and outlook
The three paleoclimate experiments, 6 ka, LGM, and LM, proposed in CMIP5/PMIP3 were performed using MIROC-ESM.Overall, MIROC-ESM gave qualitatively reasonable results for all three cases.Special care was required when conducting experiments with conditions significantly different from modern simulations, especially using a highly complex model like ESM.
The 6 ka experiment was relatively simple, requiring spinup only for the land carbon cycle.The results showed reasonable seasonal changes in the temperature distribution in response to the insolation change.The precipitation and the circulation changes suggest the enhancement of the summer African and Asian monsoons.These changes are consistent in sign with the paleoenvironmental proxy records, but the enhancement over the Sahara Desert was not large enough, as also seen in previous modelling studies (e.g.Joussaume et al., 1999;Braconnot et al., 2007).
The LGM experiment required several steps before starting the simulation.Large change in topography requires a careful adjustment to avoid the numerical instabilities by gravity waves, and very long integration for spin-up is required for the ocean physical state and carbon cycle variables.In this experiment, ocean values are taken from a LGM experiment using the previous version of MIROC AOGCM, and the spin-up for the carbon cycle was performed separately.In the results, changes in temperature and precipitation are mostly consistent with the SSTs and the pollen proxy records (MARGO Project Members, 2009;Bartlein et al., 2010).AMOC was stronger, which is contrary to the results from some reconstructions.This discrepancy is an issue for some climate models (Otto-Bliesner et al., 2007), which should be discussed in another study.
In the LM experiment, the same settings as the PI experiment were used for the fixed boundary conditions.As it is a transient experiment, preparations for time-dependent forcing data are required.Since the experiment was performed as "CO 2 prognostic", the time series of carbon emission data and initialization of the atmospheric CO 2 are also required.After the removal of the linear-regressed temperature trend, the result shows general agreement with the reconstructed hemispheric temperature variations.Since this experiment depends on solar and volcanic forcing data, updating and cross checking of the data, both for the forcing and proxy records, are continuously required.
It is again recognized that a careful and sufficiently long spin-up is mandatory to achieve an equilibrium state, especially for such a long transient simulation as LM.In ESMs, more processes in land and ocean biogeochemistry are considered, and typical timescales are different among them.In addition to the physical model part, an initialization of the carbon cycle part is also necessary.
In this paper, we have documented the implementation of the CMIP5/PMIP3 experimental protocol for paleoclimate simulations using MIROC-ESM, with some overview of the results.We have successfully performed the experiments as described here, and the results are available via the CMIP5 database (Taylor et al., 2012).Further analysis of each setting, as well as multi-model ensemble analysis, will follow this paper to contribute to the CMIP5 modelling activity.
2. As Published by Copernicus Publications on behalf of the European Geosciences Union.T. Sueyoshi et al.: Set-up of the PMIP3 paleoclimate experiments the carbon dioxide plays a crucial role in controlling radiation, ESMs are being increasingly used in global climate change projection.

Fig. 3 .
Fig. 3.The procedure of initial data preparation: (a) pre-industrial control: PI; (b) mid-Holocene: 6 ka; (c) Last Glacial Maximum (LGM): 21 ka; (d) Last Millennium: LM.Rounded rectangle denotes the change in boundary condition; shaded rectangle with arrow denotes the change in restart (initial) data, such as value overwriting and offline spin-up.For 6 ka, the integration time from the point of initialization (i.e."branch" time) was kept the same as PI experiments (dashed blue line).

Figure 4 .
Figure 4.The global annual mean sea surface temperature representation (°C) of PI in comparison with World Ocean Atlas SST (WOA, 2009).

Fig. 4 .
Fig. 4. The global annual mean sea surface temperature representation ( • C) of PI in comparison with World Ocean Atlas SST (WOA, 1998).

Fig. 6 .
Fig. 6.The change in incoming shortwave radiation (W m −2 ) at the top of the atmosphere for (a) 6 ka PI and (b) LGM-PI.

Fig. 12 .
Fig. 12.The (a) volcanic and (b) solar forcing, (c) annual mean surface air temperature in the Northern Hemisphere, (d) modelpredicted (black) and reconstructed (blue) atmospheric CO 2 concentration, and (e) detrended time series of mean annual surface air temperature anomaly (anomaly from the 1971-2000 mean).The blue shading represents the range of temperature reconstruction between the 20th and 80th percentile by Frank et al. (2010).

Table 1 .
The boundary conditions for each experiment.