The Importance of the Orbital Parameters for the Indian Summer Monsoon During the Mid-Holocene, as Deciphered From Atmospheric Model Experiments

Proxy and model-based studies suggest multi-scale temporal variability in the Indian summer monsoon (ISM). In this study, using the CESM1 atmospheric general circulation model, we carried out multiple ensemble AGCM simulations for the Mid-Holocene (MH; ≈ 6 kyr BP), Medieval Warm Period (MWP; ≈ 1 kyr BP), Little Ice Age (LIA; ≈ 0.35 kyr BP), and Historical (HS; ≈ CE 2000) periods. We used the PMIP3/CMIP5 boundary conditions for this purpose. Our simulations indicate that the ISM during the MH was stronger compared to HS and the rainfall higher, in agreement with several proxy studies. The experiments also suggest that the ISM rainfall (ISMR) was higher during MWP relative to the LIA in agreement with our earlier results from the PMIP3 models. A relatively northward migration of the ITCZ over the Indian region and strengthening of the neighboring subtropical high over the northwestern Pacific, both associated with stronger insolation associated with the obliquity and precision during the MH, seem to be important reason Indian summer monsoon during the MH.


INTRODUCTION
The Indian monsoon system is a complex phenomenon. Its seasonal evolution and variability on multiple time scales involve individual and coupled land, ocean, and atmospheric processes, in addition to the annual progression of the Earth's revolution around the Sun. From June to September, the country receives roughly 75% of the annual rainfall (Supplementary Figure 1). The seasonally phase-locked rainfall and associated circulation are referred to as the southwest monsoon or Indian summer monsoon (ISM). Any changes in the average seasonal rainfall patterns have a profound effect on agriculture, and therefore the livelihood of 1.5 billion people in South Asia (Gadgil and Rupa Kumar, 2006). The Indian summer monsoon rainfall (ISMR) is weakening over several regions (e.g., Krishnan et al., 2020), with a simultaneous increase in the extreme events in the recent ∼150 years (Boyaj et al., 2020, and the references therein). Observation-based studies in the last 50-100 years suggest that, on interannual time scales, tropical oceanic phenomena such as the ENSO dominantly influence the ISM variability (see Webster et al., 1998, and the recent reviews by Ashok et al., 2019;Mohanty et al., 2020, and the references therein). However, Proxy studies indicate that ISMR also varied on the decadal, centennial, and millennial to multi-millennial time scales, in addition to intra-annual, and interannual scales (e.g., Ramesh et al., 2010;Chakraborty et al., 2012).
The relevance of internal climate variabilities, such as ENSO or multi-decadal variations in the tropical pacific circulation, is shown to be important even in the last millennium, through analysis of PMIP3 outputs (e.g., Tejavath et al., 2020).
On the other hand, external forcings are suggested to play major roles in manifesting various past climate periods. For example, changes in the strength of monsoon over the glacialinterglacial transition are suggested to be due to changes in orbital parameters (Kutzbach and Otto-Bliesner, 1982;Kutzbach and Guetter, 1986). Earth's precession cycle of 23 kyr plays a dominant role in modulating the tropical precipitation by changing the seasonal and meridional distribution of incoming solar radiation (Kutzbach, 1981;Pokras and Mix, 1987). Speleothem reconstructions suggest that changing precession leads to antisymmetric precipitation patterns in the southern and northern hemispheres (Wang et al., 2007). Several proxy and modeling studies show that the strength of the monsoon has changed in proportion with changes in insolation on orbital scales (Kutzbach, 1981;Tuenter et al., 2005;Wang et al., 2007;Kutzbach et al., 2008;Weber and Tuenter, 2011;Shi, 2016;Kathayat et al., 2017).
The distinct external forcings have been said to have manifested the MH climate in India. A study by Gupta et al. (2005) shows that slight changes in solar radiation have brought changes in the tropical monsoon systems during the Holocene. Studies also show that Earth's climate is sensitive to small changes in solar radiation on the centennial to the millennial time scale during the Holocene (Rind and Overpeck, 1993;Shindell et al., 2001;Hu et al., 2003). Signatures of changes in insolation are reflected in the changes of the northern hemisphere monsoon systems (Staubwasser et al., 2003). Past climate records based on proxy reconstructions from India and neighborhood surroundings of the Indian region (see the review by Dixit and Tandon, 2016;Banerji et al., 2020;Tejavath et al., 2020 for the details and references) report centennial to millennial-scale changes in the ISMR during the Holocene ( 12 kyr BP to present). These studies, in general, show that from the early Holocene period (9 kyr BP) to late Holocene (3 kyr BP), the ISMR shows a long-term decreasing trend. A single model study by Kumar et al. (2019) and a multi-model study by Tejavath et al. (2020) using PMIP3 (Paleoclimate Modelling Intercomparison Project 3) show that during the mid-Holocene (MH; 6 kyr BP), the ISM was stronger relative to the present day. An analysis of the PMIP3 simulations by Tejavath et al. (2020) suggests that a robust large-scale convergence over the Indian region at 850 hPa may have been the reason for strong ISM during the MH. A study by Cobb et al. (2013) purports a highly variable ENSO during the Holocene period. However, several other proxy reconstruction studies (e.g., Gill et al., 2016Gill et al., , 2017 suggest that there was reduced ENSO activity during the MH period, and dominated by La Niña type SSTs over the tropical eastern pacific. There is another argument that the ENSO activity during the MH was indeed relatively subdued (An and Choi, 2014;Chen et al., 2019), and that such suppression is due to the enhanced Asian monsoon activity (Chen et al., 2019;Crétat et al., 2020). An ocean-atmospheric coupled model study by Crétat et al. (2020) shows orbitally driven trends during the MH to the current day (i.e., 0 kyr BP). The study also claims an increased influence of ENSO on ISMR and Indian Ocean Dipole (IOD) during the MH.
Studies based on paleo reconstruction data identify two significant periods in the last millennium (LM). These two periods, commonly known as Medieval Climate Anomaly (MWP; CE 950-1350) (e.g., Lamb, 1965)-a relatively warmer period, followed by a relatively more cooling period referred to as the Little Ice Age (LIA, CE 1500-1850) (e.g., Grove, 1988). Studies based on proxy and modeling show wetter and stronger ISM during MWP, and drier and weaker ISM during LIA relative to the LM mean (see Dixit and Tandon, 2016;Tejavath et al., 2019 for further references).
Various mechanisms have been proposed for ISM variability during the Holocene period, in addition to the changes in the orbital parameters (e.g., precession). For example, during the early Holocene climatic period, the mean latitudinal position of the Intertropical convergence zone (ITCZ) shifted northwards, resulting in an increase in the ISMR in response to the changes in the solar insolation (Fleitmann et al., 2007). Early climate model sensitivity experiments show that the timing of Northern Hemisphere (NH) summer monsoon strength has, or should have the near-zero phase relative to maxima of NH precession-driven radiation (Clemens and Prell, 2007). A proxy reconstruction study by Dixit et al. (2014a,b) shows that solar insolation played a significant role in strengthening the ISM during the early Holocene period. A study by Polanski et al. (2012) based on proxy and modeling shows that solar insolation changes strongly influenced the MH climate. Several other modeling studies suggest enhancing the northern hemisphere monsoon systems due to changes in insolation (Braconnot et al., 2007a,b;Bosmans et al., 2012;Marzin et al., 2013;Zheng et al., 2008). A recent modeling study suggests that a decrease in Saharan vegetation may have resulted in droughts in northernmost parts of India and excessive rains in southern India (Griffiths et al., 2020). A study by Leuschner and Sirocko (2003) shows that precessional and obliquity cycles coincide with the major events observed in the insolation-based ISM index. Interestingly, a study by Li and Harrison (2008) suggest in their coupled model simulations, the orbital forcing mostly enhances the Asian summer precipitation but in contrast, the neighboring ocean reduces the orbitally-induced summer precipitation, leading to the postponement of the time of summer monsoon onset over the Asian monsoon region.
However, the above studies including the coupled modeling studies, have not specifically examined the relevance of the external forcings to the ISM during MH in detail. Furthermore, the relative variations in the ISM during MH and LM have been less-studied, particularly from the modeling perspective. Keeping this in mind, using the community atmospheric model version 5 (Hersbach et al., 2015), we carry out several timeslice experiments for MH, MWP, LIA, and Historical (deemed as "present-day") periods. Importantly, we explore the sensitivity of the simulations by changing the orbital parameters; a major external forcing said to be important for the MH climate period, as discussed in detail in the next section.
This manuscript is organized as follows. In the following section, we briefly describe the model, our experimental setup, and methodology. In section "Results, " we present our results and discussion. We present our conclusions and future scope for related studies in the final section, which is section "Conclusions and Scope for Future Studies."

EXPERIMENTAL SETUP AND METHODOLOGY
The Community Atmospheric Model Version 5 (CAM5), is the atmospheric component of the well-known Community Earth System Model Version 1.2.0 (CESM1.2.0; Hersbach et al., 2015). The CESM1.2.0, formerly also known as CCSM4, is a fully coupled climate model used to generate simulations of the Earth's past, present, and future climate states. For our experiments, we configured the horizontal resolution of the CAM5 at 1.9 • latitude × 2.5 • longitude grids, with 30 vertical levels. This version of CAM is comparable and compared to the CESM simulations for the Coupled Model Intercomparison Project 5 (CMIP5) and the Paleo Model Intercomparison Project 3 (PMIP3) simulations. The PMIP3 is a collective initiative endorsed by the World Climate Research Programme (WCRP) and JSC/CLIVAR working group on coupled models and the International Geosphere and Biosphere Programme (IGBP; PAGES) .
We carried out four suites of multi-ensemble control simulations of a 30-year span for each of the MH, MWP, LIA, and Historical (Present day) periods. Each suite of the control experiments in turn contains three ensembles. Each ensemble simulation starts with January initial conditions that are different from one ensemble to another and continues for a year. The external boundary conditions used in this study are period-specific (i.e., MH, MWP, LIA, and HS), similar to those mentioned in the PMIP3 website https://pmip3.lsce.ipsl.fr/, and briefly presented in Table 1. The climatological SSTs used in this study for each period are adapted from the simulation outputs from the corresponding PMIP3 simulations that were generated with CCSM4 (Figure 1). Apart from the control simulations, we have also carried out multiple ensemble simulations for the MH period by forcing the model with different orbital forcings for different periods, e.g., present-day orbital, 8.2 kyr BP orbital, and LGM orbital forcings. The forcings are mentioned in Table 2. For comparison, the same set of initial conditions have been used in all the complementary experiments.
Further to these 30-year long control and orbital sensitivity simulations, we have also carried out suites of simulations of 1year span (simulation starts in January and ends in December) with 10-ensembles each with the same set of initial conditions across various orbital forcings. The results from the 1-year length simulations are in very good agreement with those from the 30-year long simulations, as evidenced by the evolution of the seasonal cycles (Supplementary Figure 2). Therefore, we shall not discuss the results from these runs.
Apart from the control and orbital sensitivity experiments, we have also carried out two sensitivity experiments to gauge the potential relevance of the El Niños during MH, and two more to assess the impact of the La Niñas during the same period. This involved, for example in case of the La Niña sensitivity experiments, imposing the La Niña anomalies on the climatological SSTs over the tropical pacific ocean (120 • E to 80 • W and 30 • S to 30 • N) to obtain the lower boundary SSTs representing the La Niña SSTs during the MH, and in another experiment, imposing SST anomalies throughout the tropical Indo-pacific ocean (40 • E to 80 • W and 30 • S to 30 • N). The SST anomalies imposed in these SST sensitivity experiments for the MH period are shown in Supplementary Figure 3. The SST anomalies during the La Niñas (El Niños) were obtained by compositing monthly SST anomalies of all "typical" La Niñas (El Niños). These typical events are identified as those for which the magnitude of the simulated JJAS NINO3 index is above (below) one standard deviation.
In our analysis, the climatology is taken from the entire 30year runs, because the initial conditions are taken from PMIP3 CCSM4 simulation, and they are therefore spun up well. As we are starting our simulation from January, the simulated atmosphere will adjust to the lower boundary SSTs within a couple of months. As it is, such experiments for the current day monsoonal season typically start in the month of May (e.g., Ashok et al., 2001Ashok et al., , 2004Ashok et al., , 2009Ashok et al., , 2012Guan et al., 2003).
We have calculated the monthly climatological cycles of rainfall and surface temperatures over the Indian land region bounded by 66.5 • E-101.5 • E; 6.5 • N-39.5 • N. We analyze the simulated fields of velocity potential-which represent divergence, vorticity, and moisture flux convergence. The moisture flux convergence has been computed as the sum of moisture convergence and advection. The fields of velocity potential, vorticity, and moisture convergence have been computed from the simulated outputs of horizontal winds and moisture through the application of spherical harmonics using the NCL routines. The details can be found at https://www.ncl.ucar.edu/ Applications/wind.shtml. The orbital parameters are calculated using the PMIP3 protocol https://pmip3.lsce.ipsl.fr/ and are present in Supplementary Table 1.

Simulated Mean Summer Monsoon Rainfall and Surface Temperatures Over India
Before going any further, we have validated HS simulations with the observational and PMIP3 CCSM4 historical simulations and are presented in Supplementary Figure 4. We briefly validated the seasonal cycle evolution of rainfall and surface temperatures of these simulations (PMIP3 CCSM4 and AGCM CAM5) with the ERA-20CM skin temperature (SKT; Hersbach et al., 2015), Climate Prediction Center (CPC) Global Land Surface Air Temperature data (Fan and van den Dool, 2004) and the India Meteorological Department (IMD) gridded rainfall (Rajeevan et al., 2006). We present the seasonal cycle validation in Supplementary Figure 4. From Supplementary Figure 4, we can see that the seasonal cycle evolution of rainfall and surface temperatures are well captured in the AGCM HS.

Mid-Holocene
Figures 2A,C show that simulated summer monsoon rainfall during MH is higher than the HS. Note that we compute the areaaverages over the Indian land region only (e.g., see the domain in Figure 2C). This simulated higher ISMR during MH, particularly around ∼ 6 kyr BP, agrees well with several coupled model simulations (Kumar et al., 2019;Tejavath et al., 2020), as well as with the proxy-based studies (Rawat et al., 2015;Band et al., 2018). The simulated area-averaged surface temperature found to be coldest during MH compared to the Present-day period ( Figures 2B,D), as also indicated by several coupled models of the PMIP3 vintage in Tejavath et al. (2020). The simulations suggest, during the MH, a relatively higher rainfall relative to the HS in the core monsoon region, northeast India, and foothills of the Himalayas (Figure 2C), and lower rainfall in the western part of India. Interestingly, despite being a major ISM rainfall region, the simulated rainfall over the western ghats is relatively less rainfall during the MH compared to the HS (Figure 2C).
We have observed a similar kind of pattern in the southern part of western ghats (i.e., Kerala) in PMIP3 CCSM4 coupled model simulations also. This could be because, during the MH, the monsoon winds are weaker over the western ghats region and less moisture availability compared to the HS. While the monsoon circulation seems to be stronger relative to the HS during the MH (Figure 5), when the Indus civilization, encompassing the modern-day Rajasthan and neighboring Pakistan, etc., is said to have thrived. However, as can be conjectured from the 850 hPa circulation, the low-level circulation is weaker relative to the historical period. This of course may be also model-specific. We shall explore this aspect further making a detailed PMIP model analysis and sensitivity experiments with multiple AGCMs. Apart from them, comparatively low-resolution simulations and the atmospheric-only simulation could be the reasons behind it.

MWP and LIA
Figures 3A,C, show that the area-averaged summer monsoon rainfall during MWP is higher than LIA, in agreement with a majority of the PMIP3 simulations (e.g., Tejavath et al., 2019). Figure 3B confirms that the AGCM simulations reproduce the expected surface temperature response to increased greenhouse gases, from the perspective of increased GHGs in the present day when industrialization has begun (also see Table 1), and hence considered reliable. The simulated area-averaged surface temperature during the LIA is cooler compared to the MWP, again as simulated by the coupled models of the PMIP3 vintage . Interestingly, the rainfall difference over northeast India indicates a dipolar pattern during the MWP and LIA ( Figure 3C). When area-averaged, these values are almost negligible, in agreement with results from three relatively high resolution PMIP3 simulations suggest that the area-averaged summer monsoon rainfall over northeastern India has not changed across the MWP, LIA, and historical periods (Ashok et al., 2020, under review).

Role of Tropical Ocean-Atmospheric Coupled Processes or the Lack Thereof
Apart from these, we compare the results from our simulations with the results from the corresponding past climate simulations for MH with the difference to HS. The dotted region represents a statistically significant region at a 95% confidence level from a two-tailed Student's t-test.

FIGURE 3 | Comparison between the area-averaged seasonal cycle of Medieval Warm Period (MWP) and Little Ice Age (LIA) period, (A) is for the simulated rainfall, and (B)
is for the simulated surface temperature over Indian land region. Spatial distributions of the simulated summer monsoon rainfall (mm/day; C) and surface temperature ( • C; D) for MWP with the difference to LIA. Hatched region represents a statistically significant region at an 80% confidence level from the two-tailed Student's t-test. by CCSM4 (Supplementary Figure 5), which are available as PMIP3 and CMIP5. Note that, these CCSM4, which have CAM as the atmospheric component, are fully coupled oceanatmospheric simulations (Schmidt et al., 2012;Taylor et al., 2012). From Supplementary Figure 5, we see that simulated rainfall and surface temperature climatological cycles from our experiments with the CAM AGCM are comparable to those from the CCSM4 in terms of the evolution of the annual cycle. The AGCM simulations are relatively drier and warmer over India during the summer (Supplementary Figures 5A-D) across all the climate periods compared to the coupled CCSM4 simulations. Similarly, we see a warmer bias in AGCM surface temperatures compared to the coupled model simulations. The absence of air-sea coupled processes could be a probable reason behind the dry bias.

Mid-Holocene
During the MH period, we see a stronger low pressure in the monsoon trough region northward from around 20 • N relative to the HS ( Figure 4B). This indicates the northward migration of the ITCZ, which is also seen in simulations by a few other AGCMs and coupled models of older vintage (e.g., Zhao and Harrison, 2012); strengthening of the subtropical high over the western pacific also apparently strengthens the monsoonal circulation over northern India (Figure 4B). Apart from these, the Indian subcontinent experienced relatively stronger lowlevel convergence associated with stronger monsoonal circulation and during MH ( Figure 5A). Other factors such as a stronger zonal seasonal mean SST gradient in the tropical Indian Ocean, with a warmer (cooler) eastern (western) tropical Indian Ocean by 0.1∼0.2 • C ( Figure 4A) relative to the HS, which provides a positive IOD-like background, may have also contributed to the stronger mean ISM circulation and summer monsoon rainfall. Large-scale moisture convergence is seen during the MH (Figure 5B), consequently giving a relatively enhanced summer monsoon rainfall in the core monsoon region. Western India is associated with a relatively lower moisture convergence, manifesting as lower rainfall (also see Figures 2A,C). Earlier modeling studies also show that changes in the large-scale circulation facilitates the precipitation changes during the MH period and last millennium (e.g., Polanski et al., 2014;Tejavath et al., 2019;Ashok et al., 2020). Interestingly, the simulated pre-monsoon (MAM) land temperatures during the MH are cooler than the HS (Supplementary Figure 6). A similar kind of signal can be seen in general, the PMIP3 coupled simulations also (figures not shown). From this, it seems the pre-monsoon land warming may not be relatively so important. The landsea gradient during the MH from our control run is also weak, suggesting that this may not be a major factor during the MH.

MWP and LIA
From Figure 6A, we clearly see that the Indian subcontinent experienced stronger low-level convergence, which is associated with stronger monsoonal circulation during the MWP compared to the LIA. The MWP climate period is associated with strong magnitudes of moisture convergence compared to the (C) Shows spatial distributions of simulated 500 hPa relative vorticity (×10 −5 S −1 ) for the MH period relative to that for the HS. Hatched region represents a statistically significant region at a 95% confidence level from the two-tailed Student's t-test.
LIA ( Figure 6B). These all resulted in a stronger ISM giving surplus rainfall over the Indian region during the MWP compared to the LIA (Figure 6C). (C) Shows spatial distributions of simulated 500 hPa relative vorticity (×10 −5 S −1 ) for the MWP period relative to that for the LIA. Hatched region represents a statistically significant region at a 95% confidence level from the two-tailed Student's t-test.

Possible Mechanisms for the Stronger Mid-Holocene Summer Monsoon
In Figure 7, we show the differences in the simulated insolation between the MH and the HS climate periods, associated with the changes in the obliquity. From the perspective of the Milankovitch cycles, the obliquity (axial-tilt) of the Earth during the MH is 24.1 • , and the HS is 23.4 • . In other words, during the MH, the solar activity was more in the Northern hemisphere compared to the southern hemisphere (Figure 7). This shows that during the MH climate period, the Indian subcontinent has received more solar insolation than during the Historical climate period because of the changes in the orbital parameters (Figure 7). A few studies (Bosmans et al., 2012;Zhao and Harrison, 2012) suggest that this increased-insolation likely led to deeper thermal lows on the land region and resulted in a higher land-sea thermal gradient because of the difference in thermal inertia between land and ocean. In general, this can be expected to lead to stronger monsoonal winds into the Indian sub-continent with enhanced moisture flow toward the land region, giving more precipitation (Bosmans et al., 2012;Zhao and Harrison, 2012), at least during the onset phase when the land-sea thermal contrast is important. Indeed, the seasonal atmospheric circulation, in general, can also be expected to change as a result of the change in atmospheric energy balance prompted by the strengthening of the summer by stronger insolation (Merlis et al., 2013), and also eventually influence the rainfall. From this perspective, it is interesting to see how the Indian monsoon responded to the changes in orbital parameters, particularly for the mid-Holocene time period using the AGCM.

The Mid-Holocene ISM Response to Changes in Orbital Forcings
In this section, we examine the simulated ISM response to changes in orbital forcings through an analysis of our sensitivity experiments, which were discussed in section "Experimental Setup and Methodology." We have carried out simulations for the MH climate period with different orbital forcings (e.g., present-day orbital, 8.2 kyr BP orbital, and LGM orbital forcings) are mentioned in Table 2. The motivation behind choosing the LGM and 8.2 kyr BP orbital parameters apart from HS follows: proxy-based studies suggest that during the LGM, the ISM was weaker and drier compared to the present-day (Pattanaik, 2012;Chabangborn et al., 2013). During the 8.2 kyr event, ISM had weakened abruptly (Dixit et al., 2018). Changes in orbital parameters are believed to have played a major role in weakening the ISM during the LGM (Bowen, 2009). It must be mentioned that Atlantic teleconnections also had a significant role in the weakening of ISM during the 8.2 kyr BP (Dixit et al., 2018). By applying the orbital parameters of LGM and 8.2 BP kyr for the MH gives an idea for the role of orbital parameters during MH compared to the other climatic periods, although they may possess characteristic differences in capturing the low/high rainfall patterns to these changes. The simulated MH solar irradiance over the Indian region during the summer monsoon season decreases when the orbital parameters are replaced by those observed during HS and LGM (Figures 8A,C). Furthermore, when the 8.2 kyr orbital parameters have been invoked, the simulated solar irradiance over the Indian region  is seen to be slightly more compared to the MH control simulation (Figure 8B). The other simulated changes relative to these changes in orbital parameters are discussed in the following subsections.

Simulated Rainfall and Surface Temperature
Changing the orbital forcing to that of the HS, and to that of the LGM, reduces the simulated summer monsoon rainfall relative to the MH-control experiment ( Figure 9A). Interestingly, the MH simulations with the orbital forcings pertaining to the 8.2 kyr BP resulted in relatively higher summer monsoon rainfall, contrary to proxy-based observations (Dixit et al., 2018). An earlier coupled model study also suggests that simulating the signatures of the 8.2 kyr BP event seen in proxy-observations is a difficult task (LeGrande and Schmidt, 2008). This tells us that apart from the orbital parameters, the North Atlantic teleconnections FIGURE 10 | Spatial distribution of the simulated 850 hPa differences in the time-averaged JJAS velocity potential 'χ 850' (m 2 s −1 ; Shaded), and overlaid by the mean monsoon winds (m/s; Vectors) for MH sensitivity experiments relative to that for the MH control simulation. (A) MH sensitivity experiment with HS orbital parameters to MH control simulation, (B) MH sensitivity experiment with 8.2 kyr BP orbital parameters to MH control simulation, and (C) MH sensitivity experiment with LGM orbital parameters to MH control simulation. (Dixit et al., 2018) and, importantly, meltwater forcings (Renssen et al., 2001;Wiersma and Renssen, 2006;LeGrande and Schmidt, 2008) may have played a crucial role in the weakening of ISM 8.2 kyr BP event as seen in paleo-data based studies. In this context, the simulated relatively high ISMR with the 8.2 kyr BP may be due to a lack of suitable representation of the fresh-water perturbations, as our experiments are only AGCM-based.

Simulated Circulation Changes Over India
The changing of the orbital parameters to the current day, and in a parallel experiment to that in the LGM time period, caused the weakening of the mean monsoonal winds. This resulted in weakening large-scale low-level convergence over the Indian region (Figures 10A,C), and a deficit summer monsoon rainfall, particularly over the core monsoon region (Figure 9). From Figures 9, 10, it is evident that the simulated ISMR during the MH has responded differently to changes in solar irradiance induced by changes in the orbital parameters. This suggests that enhanced solar insolation due to favorable orbital parameters during the MH time period played a crucial role in strengthening the ISM, compared to the HS and LGM, also supported by Crétat et al. (2020) This also suggests that changes in orbital parameters had a major role in the strengthening of ISM during the MH than the equatorial ocean-atmosphere dynamical processes at that time.

Potential Role of Concurrent ENSO
The simulated model response in the experiments carried out with El Niño-type of SSTs in the tropical Pacific for the MH period [Experiment ELP(MH)] resulted in negative summer rainfall anomalies over India (Supplementary Figure 7). However, the model response to the La Niña-like lower boundary forcing in the tropical Pacific [Experiment LNP(MH)] for the same period seems to be anomalous (Supplementary Figure 7), as it produces a negative summer monsoon rainfall anomaly over the Indian region. To be sure, studies such as that by Annamalai et al. (2005) suggest that some models may also need the ENSO-associated SST forcing in the tropical Indian Ocean to produce a realistic response over the Indian monsoon rainfall. In addition, recent research by Chowdary et al. (2017) suggests that the ENSO signal in the tropical Indian Ocean, in addition to its signal in the tropical Pacific may be relevant for the current day monsoon variability. In our model too, for the historical period, the Indian Ocean SST anomalies associated with La Niña, along with the Pacific SST anomalies, seem to produce above-normal summer monsoon rainfall at least from August for the historical period (Supplementary  Figure 7). Interestingly, even when we introduce the La Niñarelated SSTA in the tropical Indian Ocean to the tropical Pacific SST signals [LNPI(MH)] for the MH, the model does not simulate the above positive summer monsoon rainfall anomalies over India (Supplementary Figure 7). On the other hand, similar experiments carried out with the ENSO-SSTs in both basins simulate positive summer monsoon rainfall anomalies over India for the HS, MWP, and LIA at least from August through September.
All the exercises suggest that, at least in our model simulations, there may not be a contribution of the MH La Niñas to the wet MH summer monsoon rainfall. Their contribution to the ISM during MWP and LIA, however, is at least qualitatively palpable, in agreement with the PMIP3 results .

CONCLUSION AND SCOPE FOR FUTURE STUDIES
While several proxy studies, such as Rawat et al. (2015) and Band et al. (2018) suggest that the ISM during the mid-Holocene (MH; ∼6000 yr BP) was characterized by a wet climate relative to the modern time. The few modeling studies (Kumar et al., 2019;Tejavath et al., 2020, etc.) also suggest that the ISM was indeed stronger at that time. Interestingly, proxy-based (Mukherjee et al., 2016) and model-based studies suggest that external factors such as changes in the orbital parameters, changes in solar forcing induced by volcanic eruptions, land surface, and vegetation changes, etc., may have played an important role in evolving a distinct climate of the Earth during MH relative to the current day (Crétat et al., 2020). The present study, through various experiments with the Community Atmospheric Model (CAM5), examines the potential role that the distinct orbital parameters during the MH may have played in the manifestation of mean ISM conditions during the MH.
We have carried out 30-year long multiple ensemble simulations for MH, MWP, LIA, and Historical period (Presentday; HS) climate periods. The results show that the ISMR during the MWP is higher than that of the LIA, but lesser than that of the HS, in agreement with various PMIP3 coupled simulations . Our AGCM experiments simulate relatively higher surplus ISMR during the MH period relative to the HS, supporting proxy observation-based studies and from few modeling studies.
Results from our novel experiments show that the higher ISMR simulated in the MH is due to a stronger monsoon circulation, northward migration of the ITCZ, strengthening of the subtropical high over the western pacific relative to the HS. To examine the potential role of orbital parameters on the ISM variability during the MH, we carried out several sensitivity simulations for the MH period by changing the orbital parameter with those of the HS, LGM, and 8.2 kyr BP. These sensitivity experiments demonstrate that the higher orbital parameter during the MH, which meant more solar insolation in the northern hemisphere played a significant role in enhancing the ISM, and consequently, high ISMR. This suggests that the higher insolation during the MH, around 6 kyr BP, associated with relatively favorable external forcings (e.g., orbital parameters) likely played a significant role in enhancing the ISM during the MH compared to the HS. Of course, this is by no means a complete explanation for the high MH rainfall simulated and recorded in several proxies. For example, we have not carried out any experiments to ascertain the potential role of land-surface and vegetation changes in the mid-Holocene, which is said to have facilitated the mid-Holocene climate in West Africa (Messori et al., 2018;Crétat et al., 2020;Griffiths et al., 2020). Such experiments need a model with an interactive land-surface and vegetation model. The other factors, such as the internal variability such as ENSO, whether generated due to internal coupled dynamics or forced by orbital forcings may have played their roles. In this context, a few sensitivity simulations carried out by us suggest that our AGCM simulates below normal ISM rainfall during the MH as a response to concurrent El Niño-like SSTs imposed in the tropical Pacific (e.g., Ashok et al., 2004). This is similar to the current day ENSO impact on ISM. But we elicit a relatively low Indian summer rainfall response during the months of June and July when we force the model La Niña type of SSTs in the tropical Pacific ocean or even when the La Niña-associated concurrent SST anomalies are also imposed in the tropical Indian Ocean. Importantly, similar La Niña experiments for the MWP and LIA result in a position ISMR anomaly. This suggests that, at least in our simulations, the La Ninas have not significantly contributed to the relatively wet ISM during the MH, and that the orbital forcings may be relatively more important.
To get a comprehensive idea of the relative roles of various internal and external factors and their influences on the ISM during the MH, a more extensive study with a fully coupled landbiosphere-ocean-atmospheric model will be carried out soon. We also plan to revisit the above results by repeating these experiments with multiple AGCMs in order to ascertain our claims in this paper.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CT took lead in carrying out the simulations, analyzed the data, and wrote the manuscript. KA and SC supervised and cowrote the manuscript. All authors contributed to the article and approved the submitted version.