Seasonal modulation of the Asian summer monsoon between the Medieval Warm Period and Little Ice Age: a multi model study

Instrumental and proxy records indicate remarkable global climate variability over the last millennium, influenced by solar irradiance, Earth’s orbital parameters, volcanic eruptions and human activities. Numerical model simulations and proxy data suggest an enhanced Asian summer monsoon during the Medieval Warm Period (MWP) compared to the Little Ice Age (LIA). Using multiple climate model simulations, we show that anomalous seasonal insolation over the Northern Hemisphere due to a long cycle of orbital parameters results in a modulation of the Asian summer monsoon transition between the MWP and LIA. Ten climate model simulations prescribing historical radiative forcing that includes orbital parameters consistently reproduce an enhanced MWP Asian monsoon in late summer and a weakened monsoon in early summer. Weakened, then enhanced Northern Hemisphere insolation before and after June leads to a seasonally asymmetric temperature response over the Eurasian continent, resulting in a seasonal reversal of the signs of MWP–LIA anomalies in land–sea thermal contrast, atmospheric circulation, and rainfall from early to late summer. This seasonal asymmetry in monsoon response is consistently found among the different climate models and is reproduced by an idealized model simulation forced solely by orbital parameters. The results of this study indicate that slow variation in the Earth’s orbital parameters contributes to centennial variability in the Asian monsoon transition.


Introduction
The last millennium, from 850 A.D. to the present, is a key period in understanding Earth's climate response to external forcing and its internal variability on decadal-to-centennial timescales (e.g., Crowley 2000;Jones et al. 2009;PAGES 2k Consortium 2013;Smerdon and Pollack 2016). A variety of instrumental and proxy records indicate that the global climate has experienced remarkable fluctuations during this period due to large volcanic eruptions (e.g., Miller et al. 2012;Schurer et al. 2014;Yang and Jiang 2015;Liu et al. 2016), variations in total solar irradiance (TSI) (Gray et al. 2010;Lean 2010), Earth's orbital parameters (Berger and Loutre 1991;Kaufman et al. 2009), greenhouse gases (GHGs) (e.g., Schmidt et al. 2011), and internal climate variability. On a centennial timescale, anthropogenic global warming accompanying rapid GHG increases since the 19th century is the most remarkable climate fluctuation (e.g., Bradley et al. 2003;Jones and Mann 2004;Ammann et al. 2007;Phipps et al. 2013;Masson-Delmotte et al. 2014). In addition, predominant global climate perturbation from the Medieval Warm Period (MWP) of 950-1250 A.D. to the Little Ice Age (LIA) of 1400-1700 A.D. (e.g., Mann et al. 2009) has also attracted much attention because it is the most recent example of long-term, naturallygenerated global climate variability (e.g., Liu et al. 2011;Man et al. 2012;Otto-Bliesner et al. 2016;Shi et al., 2016a, b).
Historical variations in the Asian monsoon have important links to environmental changes and the cultural development of human societies over South and East Asia (Clift and Plumb 2008;Fu et al. 2008). Past monsoon variability has been addressed both from paleoenvironmental reconstructions based on proxy records and from numerical climate simulations (e.g., Kutzbach et al. 2008;Zhang et al. 2008;Battisti et al. 2014;Masson-Delmotte et al. 2014;Fallah and Cubasch 2015;Yan et al. 2015). Variations in summertime temperature over the continents in the Northern Hemisphere (NH) in response to seasonal insolation anomalies brought on by changes in orbital parameters (including the longitude of perihelion precession) are an essential factor in past monsoon variability (e.g., COHMAP members 1988; Mohtadi et al. 2016). The land-sea contrast in tropospheric temperature favored by summertime insolation is further amplified by latent heating of the monsoon rainfall that establishes deep-tropospheric monsoon circulation (e.g., Li and Yanai 1996;Webster et al. 1998;Dai et al. 2013;Wu et al. 2015;Zhang et al. 2016). On a millennial timescale, past monsoon variability has largely been controlled by precession variability (e.g., Kutzbach et al. 2008;Merlis et al. 2013;Mohtadi et al. 2016) and amplified by obliquity variability (Bosmans et al. 2015;Wu et al. 2016) as it changed the large-scale thermal contrast between land and sea. Oxygen isotope-based reconstructions of Indian and East Asian monsoon intensity show coherent variability on a millennial timescale (e.g., Battisti et al. 2014;Kathayat et al. 2016;Mohtadi et al. 2016), suggesting that orbital forcing has a strong effect on continental-scale monsoon variability.
On a centennial timescale, orbital forcing also played an important role in long-term differences of the intertropical convergence zone and associated monsoon activity between the MWP and LIA (e.g., Sachs et al. 2009;Oppo et al. 2009;Schneider et al. 2014;Yan et al. 2015). Major characteristics of the monsoon during these two periods have been reconstructed by proxy records that are broadly available across Asia, especially in China (e.g., Chen et al. 2015). Historical documents and speleothem records suggest that precipitation over northern China associated with the East Asian summer monsoon increased significantly during the MWP and decreased during the LIA (Ren and Zhang 1996;Wang et al. 2001;Qian et al. 2003;Tan 2007). Oxygen isotope data retrieved from Indian stalagmite records also suggest an enhanced Indian summer monsoon during the MWP (Sinha et al. 2011(Sinha et al. , 2015Berkelhammer et al. 2014). The coherent enhancement of the East Asian and Indian summer monsoon during the MWP, as compared to the LIA suggests that orbitally forced summertime insolation variability is also a key factor for centennial-scale monsoon variability during the last millennium.
Monsoon variability over the last millennium has also been addressed by numerical simulations based on single (Fan et al., 2009b;Liu et al. 2009Liu et al. , 2011Man et al. 2012;Polanski et al. 2014) or multiple climate models (Fallah and Cubasch 2015;Rojas et al. 2016;Shi et al., 2016a, b;Yang and Jiang 2017) estimated as part of the Paleoclimate Modeling Intercomparison Project Phase 3 (PMIP3) ). The climate models reproduced, in large part, the enhanced Asian summer monsoon and increased precipitation over northern China (35°N-42°N, 100°E-120°E) during the MWP (e.g., Man et al. 2012;Shi et al., 2016b), consistent with the proxy evidence. The enhanced Asian summer monsoon during the MWP, as compared to the LIA, is associated with larger land-sea temperature contrasts (e.g., Man et al. 2012;Shi et al., 2016a, b) induced by strong summertime insolation ( Fig. 1; detailed in the Methods section). These large-scale features of the Asian summer monsoon during the MWP are consistent among multiple climate model simulations (Shi et al., 2016a) with nonnegligible inter-model spread in regional rainfall patterns (Shi et al., 2016b).
In addition to seasonal-mean features, seasonal transitions (changes in onset and withdrawal timing) are also important aspects of paleo-monsoon variability (e.g., Fig. 1 Seasonal variation of insolation (W m −2 ) at the top of the atmosphere averaged over 10 PMIP3 models. Shading represents the anomaly between the MWP (950-1250 A.D.) and LIA (1400-1800 A.D.; hereafter ML anomaly). Contours represent climatology in the LIA Ueda et al. 2011;Shi 2016;Kong et al. 2017). The peak NH insolation anomaly due to the Earth's orbit variability is not necessarily concurrent with the summer solstice. As for the difference between the MWP and LIA, time-varying precession (detailed in the Methods section) leads to a positive NH insolation anomaly with its peak in August to September as shown in Fig. 1 (Kaufman et al. 2009). In contrast, a negative insolation anomaly is found in May to June, around the timing of the monsoon onset, implying a seasonally asymmetric response of the Asian summer monsoon to orbital forcing (e.g., Shi 2016; Kong et al. 2017) during the last millennium.
This study reveals that the slow variation in orbital parameters from the MWP to the LIA results in an inverse response of the Asian monsoon between early and late summer, as shown by simulations of the last millennium using multiple climate models and sensitivity experiments by an atmosphere-ocean coupled general circulation model (AOGCM). The results indicate that centennial Asian monsoon variability shown in previous works is not seasonally uniform, but rather characterized by changes in the seasonal transition of the monsoon between the MWP and LIA, which is suggestive for past, present and future climate variability.

PMIP3 multiple model ensemble
We analyzed climate simulations of the last millennium (850-1850 A.D.) of 10 AOGCMs (hereafter denoted as past1000) from PMIP3 ) and the Coupled Model Intercomparison Project phase 5 (CMIP5) (Taylor et al. 2012). Table 1 summarizes the models used in this study. Although simulation outputs from 11 models are available from the CMIP5 Earth System Grid Federation data server as of 01 February 2017, we excluded output from the MIROC-ESM due to climate drift in its simulation (Sueyoshi et al. 2013).
External forcing factors for the past1000 runs include time-varying TSI, volcanic eruptions, GHGs, land use, and orbital parameters. Details of the experimental setup are found in Schmidt et al. (2011Schmidt et al. ( , 2012. Note that the models employed different TSI and volcanic forcing data as listed in Table 1 (see Schmidt et al. 2011Schmidt et al. , 2012. This difference contributes partly to inter-model spread among simulation outputs. Varying land uses over time such as gradually increasing cropland and pasture area (Schmidt et al. 2012) prescribed in these simulations can affect regional climate variability, including surface temperature and monsoon systems (Takata et al. 2009;He et al. 2014). Historical GHG concentrations (CO 2 , CH 4 and N 2 O) are derived from ice cores in Antarctica (Schmidt et al. 2011). Orbital parameters that vary over time, including eccentricity, obliquity and precession result in low-frequency variability in the seasonal-latitudinal distribution of insolation as shown in Fig. 1 (Kaufman et al. 2009;Schmidt et al. 2011). Table 2 summarizes the orbital parameters averaged during the MWP and LIA. The negative early summer insolation anomaly in the NH during the MWP relative to the LIA is largely due to a shift of the perihelion (Schmidt et al. 2011). In this study, we compare long-term mean climates of the MWP and LIA. Note that we determined the MWP climatology in FGOALS-gl by averaging the period from 1000 to 1250 A.D. because only outputs after 1000 A.D. are provided by this model. Figure 1 shows the 10-model mean insolation anomaly (MWP minus LIA; hereafter ML) at the top of the atmosphere. The seasonal NH insolation anomaly due to orbital forcing, especially precession, is remarkable (−1.5-3 W m −2 ) despite its slow variation on a centennial timescale (e.g., the difference in perihelion is 7.28°; Table 2). A large positive anomaly in August and September (1.5-2.5 W m −2 ), in contrast to a small anomaly in June ( Fig. 1), enhances the land-sea thermal gradient and the Asian summer monsoon (Man et al. 2012). In contrast, a broad negative insolation anomaly is found in the NH in the winter-to-spring period (peaking at −1.5-−0.5 W m −2 in March to April), resulting in the strong seasonal reversal of the NH insolation anomaly. This seasonally asymmetric forcing potentially Table 1 List of PMIP3 models used in this study. The TSI and volcanic columns indicate forcing data sets prescribed in the models.   affects the seasonal transition of the Asian summer monsoon (e.g., Kong et al. 2017).

Idealized orbital forcing experiments in an AOGCM
The Asian monsoon simulated in the past 1000 runs is influenced by the time-varying external forcing listed above. The seasonally asymmetric orbital forcing ( Fig. 1) may have an especially strong effect on the Asian monsoon (Ueda et al. 2011;Man et al. 2012;Kong et al. 2017). To evaluate the contribution of the orbital forcing to the ML anomaly in the Asian monsoon, we conducted idealized simulations using an AOGCM, the Meteorological Research Institute (MRI) Coupled General Circulation Model version 2.3 (MRI-CGCM2.3) (Yukimoto et al. 2006a). MRI-CGCM2.3 has widely been used for future climate projections (Yukimoto et al. 2006b) and for paleoclimate simulations including the mid-Holocene, 6 ka Braconnot et al. 2007 (Inoue and Ueda 2009). Its atmospheric model is a spectral model with T42 spatial resolution (approximately 2.8°) that uses 30 vertical layers. The ocean component is a Bryan-Cox-type global ocean model with a horizontal resolution of 2.5°longitude and 2.0°-0.5°latitude and 23 layers. We integrated the model without any flux adjustments. The CO 2 concentration and solar constant were set to 280 ppmv and 1365 W m −2 , respectively. We set the other boundary conditions as presented in Table 1 of Kamae et al. (2016) for their pre-industrial control simulations, and integrated the model for 500 years to spin up the model in the same way as Kamae et al. (2016) had done for their pre-industrial control run. Next, we changed the orbital parameters as shown in Table 2 and integrated the model for another 55 years for the MWP and LIA cases. In this paper, we call these two simulations "orbital forcing runs." We performed an ensemble of five member simulations for each with slightly different initial conditions derived from the pre-industrial run with 5-year intervals. We determined the climatology as an average over the last 50 years for the individual runs (250 years in total). Note that this study only addresses the transient response of the Asian monsoon to radiative forcing on a multi-decade timescale. Effects of long-term deep-ocean temperature variations and overturning circulations on the Asian monsoon during the MWP and LIA are not examined in this study.

Results and Discussion
Seasonally asymmetric response in PMIP3 model simulations In this section, we compare the seasonal evolutions of the Asian monsoons during the MWP and LIA simulated in multiple PMIP3 models. Figures 2 and 3 show the ML rainfall and atmospheric circulation anomalies over the Asian monsoon region averaged for April-May-June (AMJ) and August-September-October (ASO). For spatial patterns, we chose these months to compare general features of late spring to early summer (around the monsoon onset) period, with its negative insolation anomaly (Fig. 1) and late summer to early fall (around monsoon withdrawal), with its positive insolation anomaly, and to emphasize the seasonal asymmetry in the atmospheric circulation and rainfall responses. Note that the timing of seasonal monsoon transitions varies by region (e.g., onset  Table 3 over the South China Sea in mid-May is 1 month earlier than that over Korea and Japan in early to mid-June; e.g., Tao and Chen 1987;Wang and Lin, 2002;Ueda et al. 2009). In this study, our main focus is continental-scale Asian monsoon variability including the South and East Asian monsoons. We also show seasonal evolution of the ML anomalies in individual months (detailed below).
In the 10-model ensemble mean, ASO rainfall is more active over Asia during the MWP than during the LIA (Fig. 2b), similar to June-July-August (JJA) mean anomalies shown in previous studies (Man et al. 2012;Shi et al., 2016a). The positive ML anomaly, which has high inter-model correspondence, can be found broadly over the Indian subcontinent, northeastern India, Southwest China and Eastern China. The mid-to-upper tropospheric thickness (200-500 hPa; Fig. 3b), an index for mid-to-upper tropospheric temperature (e.g., Dai et al. 2013;Kamae et al. 2017; see the Introduction section), shows general warming over the Eurasian continent, especially over the middle latitudes. In the PMIP3 multimodel mean, the summertime climatological position of the upper-tropospheric South Asian High (Shi et al., 2016a) and its seasonal evolution (intensification and northward shift from spring to summer, followed by zonal expansion in summer; e.g., Yanai and Wu 2006) are well reproduced (figure not shown). The spatial pattern of vertical wind shear (850−200 hPa) indicates enhanced low-level (upper-level) monsoon westerly (easterly) over the Arabian Sea, Indian subcontinent, Bay of Bengal, Southeast Asia, South China Sea, and Eastern China (Fig. 3b). Tropospheric temperature and associated vertical wind shear are widely used as indices for Asian monsoon intensity (Webster and Yang 1992;Li and Yanai 1996;Kawamura 1998;Kamae et al. 2014).
These anomalies correspond well with the climatological structure of the Asian monsoon ( Fig. 6a in Kamae et al. 2017), indicating a general enhancement of the Asian summer monsoon. The spatial patterns of tropospheric warming and associated anomalous atmospheric circulation and rainfall are partly dependent on which model is used (Shi et al., 2016a, b) especially over the western North Pacific (detailed below), but mid-latitude Eurasian continental warming and the enhanced monsoon are consistent among the different climate model simulations as shown in Figs. 2b and 3b (Man et al. 2012;Shi et al., 2016a, b).
In contrast to the enhanced late summer monsoon, AMJ-mean precipitation (Fig. 2a) shows a negative anomaly over Southeast Asia and Southern China and a positive anomaly over the western tropical Pacific. Figure 4a shows the climatological seasonal variation of rainfall averaged over 15°N-35°N and its ML anomaly. The model captures an earlier climatological monsoon onset over Southeast Asia (90°E-120°E in May; Tao and Chen 1987;Wang and Lin, 2002) than over the Arabian Sea (60°E-80°E in June to July). The negative AMJ rainfall anomaly over Southeast Asia is distinct from the positive ASO anomaly (70°E-120°E; Figs. 2b, 4a). Table 3 summarizes ML anomalies in seasonal rainfall increment (ASO minus AMJ) averaged over the Asian region shown by the rectangles in Figs. 2 and 4a (15°N-35°N,  70°E-120°E). Although inter-model spread of the JJAmean ML anomaly in the Asian monsoon rainfall is remarkable (Shi et al., 2016b), the seasonal asymmetry in the ML rainfall anomaly is robust among the PMIP3 models (0.09 ± 0.03 mm day −1 ). Note that an East Asian-Western North Pacific monsoon index, determined by a meridional gradient of low-level zonal wind (20°N-30°N, 110°E-140°E minus 5°N-15°N, 90°E-130°E) Page 5 of 13 , shows a large intermodel spread (0.05 ± 0.28 m s −1 ). The rainfall anomaly over East Asia is also highly dependent on which model is assumed (Shi et al., 2016b). The model diversity may be related to a large intermodel spread in the ML anomaly in the Pacific Ocean (Shi et al., 2016b). The possible importance of the Pacific Ocean to the East Asian monsoon anomaly is addressed in the Discussion section. Seasonally asymmetric responses can also be found in monsoon circulation and continental-scale temperature patterns. Figure 3a shows AMJ-mean tropospheric thickness and vertical wind shear. In contrast to ASO (Fig. 3b), a negative anomaly with its peak at 50°N is broadly found over the Eurasian continent. An anomalous meridional temperature gradient associated with this cooling results in a low-level easterly wind anomaly over East and South Asia (5°N-45°N). This easterly wind anomaly is indicative of a delay in the seasonal migration of the Asian summer monsoon (reversal of wind direction). Figure 4b shows the meridional tropospheric temperature gradient (25°N-60°N minus EQ-15°N) over Asia. Except for the AMJ period, the NH mid-latitude tropospheric temperature is higher during the MWP than during the LIA. Figure 5a shows the simulated annual-mean ML anomaly in sea surface temperature (SST). In the PMIP3 ensemble mean, the surface and troposphere over the NH middle latitudes are generally warmer during the MWP than in the LIA in response to external forcing including GHGs and orbital forcing (Schmidt et al. 2011(Schmidt et al. , 2012. Despite the generally warmer conditions, a widespread cool anomaly is found in the AMJ period (Figs. 3a, 4b), consistent with the anomalies in circulation and precipitation over Asia having opposite signs to those during ASO (Figs. 2, 3, 4). This continentalscale structure of the tropospheric temperature anomaly is a robust feature among the PMIP3 models (Table 3). Note that the East Asian-Western North Pacific monsoon index ) is largely consistent with the South Asian monsoon (stronger monsoons in ASO than in AMJ) but shows a large inter-model spread (Table 3). More details on the inter-model divergence are presented in the Discussion section.

Effects of orbital forcing
In the previous section, we showed that there is consistent seasonal reversal of sign in ML anomalies in precipitation, atmospheric circulation, and continental-scale tropospheric temperature among the PMIP3 models. Table 3 Seasonal differences (ASO minus AMJ) and 95% confidence intervals for precipitation (mm day −1 ), tropospheric thickness (between 200 and 500 hPa levels; m) over Asia, and an East Asian-Western North Pacific monsoon index, the meridional gradient of zonal wind at the 850 hPa level ) during the MWP relative to the LIA. PMIP3 (orbital forcing run) row represents a 10-model ensemble mean (5-member ensemble mean for the MRI-CGCM2.3 orbital forcing run) PMIP3 Orbital forcing run These seasonally asymmetric responses are difficult to understand as a forced response to seasonally uniform forcing that includes TSI and GHGs (e.g., Ueda et al. 2006;Endo and Kitoh 2014). Here, the seasonal NH insolation anomaly due to changes in orbital parameters (Table 2) has opposite signs in AMJ and ASO (Fig. 1), consistent with the Asian monsoon response. In this section, results of the orbital forcing runs conducted in MRI-CGCM2.3 are compared with the PMIP3 ensemble to examine the physical relationship between orbital parameters and the Asian monsoon. Figure 6 shows the results of the five-member ensemble simulations with prescribed MWP and LIA orbital parameters. The ML rainfall anomaly is largely negative over Southeast Asia and Southern China (15°N-35°N, 90°E-120°E) in AMJ (Fig. 6a). The AMJ negative anomaly is indicative of a delay in the seasonal migration of the Asian summer monsoon, similar to the PMIP3 results (Fig. 2a). In ASO, positive anomalies are found over the Indian subcontinent, Southeast Asia, and southern China (Fig. 6b). However, negative anomalies are also found over the Bay of Bengal and the Yellow River Basin, distinct from the PMIP3 ensemble mean (detailed below). The seasonal sign reversal (from negative to positive) of the rainfall anomaly averaged over the Asian region (Fig. 6) is qualitatively consistent with the PMIP3 model ensemble (Table 3). Figures 7 and 8 show the seasonal responses of tropospheric thickness and rainfall. An negative tropospheric temperature anomaly in AMJ and a positive one in ASO are broadly found over the mid-latitude Eurasian continent, similar to the PMIP3 outcomes. Associated easterly anomalies over northeastern, eastern and southern China and southern India are found in AMJ, in contrast to an ASO westerly anomaly over Japan, northeastern and eastern China, and western India. The signs of the rainfall and tropospheric temperature anomalies over Asia reverse by season from negative to positive (Fig. 8a, b), largely similar to the PMIP3 results (Figs. 3, 4). This suggests that orbital forcing is an important factor in the differing monsoon responses of the MWP and LIA.
The simulated responses to orbital forcing are somewhat different from the PMIP3 multimodel mean. The Asian rainfall anomaly (negative during AMJ and positive during ASO; Figs. 6, 8a) is generally weaker and less spatially uniform than that estimated by the PMIP3 (Figs. 2, 4a), contributing to less statistical significance ( Table 3). The simulated ML anomaly in the tropical SST is generally negative, while the PMIP3 multimodel mean shows a positive anomaly (Fig. 5). Tropospheric  Fig. 2a and b Fig. 2b, but for MRI-CGCM2.3 orbital forcing runs Kamae et al. Progress in Earth and Planetary Science (2017) 4:22 temperature has a larger negative anomaly in AMJ and a smaller positive ASO anomaly, and is less zonally uniform (Figs. 6, 7, and 8b) than in the PMIP3 (Figs. 3, 4b).
In the next section, we discuss the limitations of the sensitivity experiments we conducted to quantify effects of orbital forcing on Asian monsoon variability over the last millennium.

Possible importance of other external forcing factors
The idealized orbital forcing run examined in the previous section reproduces seasonal asymmetry in the Asian monsoon rainfall and circulation anomalies similar to the PMIP3 model ensemble mean. However, these two (the PMIP3 and orbital forcing run) also show distinct regional responses. The PMIP3 model ensemble mean produces larger ML temperature and rainfall anomalies (warmer and wetter during the MWP than the LIA in both AMJ and ASO; Figs. 2, 3, 4, 5a) than the orbital forcing run (Figs. 5b,6,7,8). Especially, the ML tropical tropospheric temperature anomaly is broadly negative (EQ-20°N) in the orbital forcing run (Fig. 7), resulting in a low-latitude circulation anomaly (e.g. westerly over northern India in AMJ) that is inconsistent with the PMIP3 (Fig. 3). The Asian rainfall anomaly is less zonally uniform than the PMIP3 estimates (Figs. 2, 6). These differences may be ascribed to less ML difference in the global-mean temperature between the MWP and LIA in the orbital forcing run than in the PMIP3 case.  also showed that the ML warming anomaly in the last-millennium orbital forcing run is smaller than in all forcing runs. The PMIP3 multimodel mean produces a larger global SST warming anomaly during the MWP than in the LIA, particularly over the middle and high latitudes (Fig. 5a), in contrast to the tropical-wide negative SST anomaly that the orbital forcing run estimates (Fig. 5b). This SST difference may make a substantial contribution to the difference in the mid-to-upper tropospheric temperatures, where we see a cooler troposphere in the orbital forcing run than the PMIP3 as in Figs. 3 and 7 (Ueda et al. 2006;Kamae et al. 2015). In   Fig. 4a and b Fig. 4b, but for orbital forcing runs Kamae et al. Progress in Earth and Planetary Science (2017) 4:22 Page 8 of 13 contrast to increased summer insolation (Fig. 1), winter insolation inversely decreases over the NH and SH, resulting in limited global-mean annual-mean insolation anomalies. This limited global-mean forcing cannot simply explain the globally-warmer climate during the MWP than in the LIA (Fig. 5a). Under the PMIP3 framework, GHGs, volcanic aerosols, TSI, and land cover all varied over time in the climate model simulations in addition to the orbital forcing. For example, the CO 2 concentration during the MWP is higher than in the LIA, especially after 1600 A.D. (Schmidt et al. 2011), contributing partly to global warming and the larger summertime tropospheric warming over the Eurasian continent than over the surrounding oceans as shown in Fig. 3b (Kamae et al. 2014). In addition, TSI and large volcanic eruptions are also important factors in the ML climate anomaly (e.g., Schurer et al. 2014;Otto-Bliesner et al. 2016). Phipps et al. (2013) showed that the lastmillennium temperature variation, in response to orbital forcing, is smaller than that in response to all forcing factors. Orbital forcing is substantially important for the seasonally asymmetric Asian monsoon response (Fig. 8) although it is a minor factor in the global ML climate anomaly (Fig. 5). The effect of volcanic eruptions on the SST anomaly is also discussed in the next subsection.

Role of sea surface temperature pattern
The PMIP3 multi-model mean reproduces general sea surface warming during the MWP than during the LIA (Fig. 5a). We should note that the spatial pattern of the ML SST anomaly and associated regional anomalies in the Asian and Western North Pacific monsoon (rainfall and atmospheric circulation) are highly dependent on the models used, as can be seen in Table 3 as per Shi et al. (2016b), who compared ML SST anomalies simulated in the individual PMIP3 models and found that (1) generally, mid-to-high latitude warming is larger than in the tropics; and (2) inter-model spread in the regional SST anomaly is substantially large over the equatorial Pacific and the mid-latitude North Pacific. A majority of the PMIP3 models simulate global ocean warming with limited east-west gradients in the equatorial Pacific Shi et al., 2016b). In contrast, one PMIP3 model, MRI-CGCM3, reproduces a La Niñalike SST anomaly, largely consistent with proxy evidence Shi et al., 2016b). This equatorial Pacific SST anomaly contributes to the simulated East Asian rainfall anomaly (Shi et al., 2016b). These results indicate that the reproducibility of the SST pattern over the tropical Pacific and the western North Pacific is essential in simulations of the Asian monsoon for the last millennium. It is also noteworthy that large uncertainties still exist in the centennial variability of the equatorial Pacific SST associated with ENSO activity and tropical volcanic eruptions during the last millennium (e.g., Adams et al. 2003). Mann et al. (2009) revealed that positive tropical radiative forcing due to high solar irradiance and inactive volcanism are responsible for the more La Niña-like mean state during the MWP as compared to the LIA, inconsistent with a more El Niño-like condition suggested in Yan et al. (2011). In the near future, we plan to conduct a set of atmospheric simulations prescribing tropical SST patterns (i.e., La Niña-like and El Niño-like patterns) to evaluate the relative importance of the SST patterns to regional and seasonal Asian rainfall and the effect of orbital forcing-induced landsea warming contrast.

Data-model and model-model comparisons
A variety of proxy records for the last millennium Asian monsoon could improve our understanding of the long-term variation in the Asian monsoon's seasonal transition. However, it is not easy to procure such records because temporally high-resolution (i.e., sub-seasonal) information is needed to evaluate datamodel consistency in the seasonal monsoon transition. One of promising data source for this purpose is meteorological information in historical documents including diaries. Adamson and Nash (2014) reconstructed regional monsoon behavior over western India for the period 1781-1860 from diaries and other historical documents. They converted qualitative information in the historical documents into a quantitative monthly rainfall index. Wang and Zhang (1992) reviewed several studies trying to reconstruct the historical variation of monthly rainfall and onset/withdrawal dates of the Meiyu rainy season in the Lower Yangtze region in the 18th century. Daily meteorological information recorded in historical documents is beneficial for reconstructing the variability of the Asian summer monsoon. However, such meteorological records were rarely continued for even several decades, let alone centuries, so they are not broadly effective in evaluating long-term monsoon variability. Further efforts to reconstruct past climate variability over Asia by combining a variety of proxy data including tree-ring records (e.g., Cook et al. 2013;Fan et al., 2009a;Fritts 2001;Takahashi et al. 2011), lake-freezing records (e.g., Arakawa 1954, Sharma et al. 2016, and phenological records (e.g. Aono and Omoto 1994;Mikami 2008) may improve our understanding of the long-term variability of the Asian monsoon.
In this study, the effect of orbital forcing is evaluated using only a single climate model. The derived results suffer from model-dependency problems including the effects of model biases. More systematic efforts are needed to better evaluate each particular forcing effect and the results' dependency on models. Previous studies have only evaluated the effects of respective forcing on continental-scale, NH-mean, or global-mean climate features (e.g., Phipps et al. 2013;Schurer et al. 2013Schurer et al. , 2014Otto-Bliesner et al. 2016). Monsoon responses to historical land use changes (e.g., Takata et al. 2009), aerosol forcing (e.g., Bollasina et al. 2011;Polson et al. 2014), and TSI (e.g., Shi et al. 2014) should be further evaluated in a multi-model framework. Further modeling work may contribute to better understanding of the forced monsoon response and its internal variability on interannual, decadal, centennial, and millennial timescales.

Conclusions
We examined long-term variation in the seasonal transition of the Asian monsoon during the last millennium using ten climate model simulations. The models reproduce an enhanced late summer MWP Asian monsoon but weakened early summer Asian rainfall compared to the LIA. The cooling-to-warming transition over the Eurasian continent before and after June with limited variation over the oceans results in the seasonally asymmetric deep-tropospheric circulation and rainfall responses over Asia. The robust seasonal asymmetry in the Asian monsoon response among the climate models is qualitatively reproduced by an AOGCM simulation forced by orbital parameters, suggesting that the seasonal-latitudinal distribution of the insolation anomaly is important to the centennial variability of the seasonal transition of the Asian monsoon.
Orbital forcing cannot explain the annual-mean climate anomaly between the MWP and LIA, implying that other forcing factors, including GHGs, have large contributions. Tropical SST response to orbital forcing in particular differs from the multiple model simulations. This inconsistency contributes to the tropospheric temperature and circulation anomaly over Asia. The east-west SST pattern in the Pacific shows remarkable inter-model variance, contributing substantially to the mid-latitude terrestrial climate anomaly between the MWP and LIA (Shi et al., 2016b). The La Niña-like ML SST anomaly in the Pacific suggested by proxy records was not successfully reproduced in the majority of climate model simulations. Uncertainty in centennial tropical SST variability and its contribution to the monsoon during the last millennium still remains a challenge. Further efforts in data-model synthesis including comparisons of simulation outputs with paleoenvironment reconstructions and paleoclimate data assimilation (e.g., von Storch et al. 2000;Fang and Li 2016) are needed to advance our understanding of the last millennium's monsoon variability.  Table 1 for producing and making their model output available. For CMIP5/PMIP3, the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provided coordinating support, and led the development of the software infrastructure in partnership with the Global Organization for Earth System Science Portals.

Funding
This work was supported by JSPS KAKENHI Grant Numbers 17 K14388 and 17 K01223.

Availability of data and materials
The PMIP3 experimental data are available from the CMIP5 Earth System Grid Federation data portal. The boundary forcing data for the PMIP3 experiments are available as supplemental data of Schmidt et al. (2011Schmidt et al. ( , 2012. The model outputs from the orbital-forcing sensitivity experiments are accessible via the authors.
Authors' contributions YK drafted the manuscript. YK, MO, and HU proposed the topic and conceived and designed the study. TK and MO analyzed the multi-model data and helped in interpretation. YK conducted the sensitivity experiments and analyzed the simulated results. All authors read and approved the final manuscript.
Authors' information MO's current affiliation is Weathernews Inc.