1 Introduction

Monsoon is a forced response of the coupled atmosphere–ocean–land system to annual variation of insolation. It is characterized by an annual reversal of prevailing surface winds and a contrast between wet summer and dry winter. The present analysis focuses on monsoon precipitation as it is a more important variable than the winds to affect the environment and society.

Based on monsoon rainfall characteristics, monsoon domains can be delineated by the regions where the annual range (local summer mean minus winter mean) of precipitation rate exceeds a threshold of 2.0 mm/day and the local summer precipitation exceeds 55 % of the annual total (Wang and Ding 2006; Liu et al. 2009). Here summer is defined as May through September (MJJAS) for the NH and November through next March (NDJFM) for the SH. The defined monsoon precipitation domain is shown in Fig. 2.

The climate variability of the global monsoon can be measured by the global summer monsoon precipitation (GSMP), which is defined as the local summer precipitation in the global monsoon domain, i.e., MJJAS precipitation in the NH and the ensuing NDJFM precipitation in the SH. Since the annual range in monsoon regions is dominated by the local summer precipitation, the GSMP represents approximately the amplitude of the annual variation of precipitation.

Future climate change involves both internal and forced variability in the climate system. To distinguish contributions to climate change from internal variability and forced variability, it requires investigation of climate variations on a longer period, such as in the last millennium. Unfortunately, such an analysis is impeded by substantial limitations of the availability and quality of global historical records.

An alternative and effective approach is to probe causes of the climate variability simulated by coupled climate models that can reproduce reliable long-term climate variations. The underlying assumption is that the models’ physics is robust enough for providing adequate physical understanding of the causes of the past climate variations. This type of millennial simulations has been constructed using complex coupled models, including the HadCM2 model (Johns et al. 1997), CSIRO model (Vimont et al. 2002), ECHO-G model (Rodgers et al. 2004), the NCAR-CSM (Mann et al. 2005; Ammann et al. 2007) and others. Various aspects of the climate variability simulated by the ECHO-G model have been examined, including the intrinsic internal variability, temperature variation, ENSO, AMOC (Atlantic meridional overturning circulation), NAO (North Atlantic Oscillation) and monsoons (Zorita et al. 2003, 2005; Gonzalez-Rouco et al. 2003; Gouirand et al. 2007a, b; Rodgers et al. 2004; Min et al. 2005a, b; Wagner et al. 2005; Liu et al. 2009). The results gained from these studies have generally built on the model’s credibility for comprehension of pertinent physical processes responsible for the simulated climate variability.

Thus far, no analysis of the secular changes of GSMP has been carried out in terms of any climate model’s millennial integration except a recent work by Liu et al. (2009). They found that the strength of the global-average summer monsoon precipitation in the forced run exhibits a minimum during the Little Ice Age (LIA) and a maximum during Medieval Warm Period (MWP), which follows the natural variations in the total amount of effective solar radiative forcing. The notable strengthening of the global summer monsoon rainfall in the last 30 years of the forced simulation (1961–1990) is unprecedented and owed in part to the increase of atmospheric CO2 concentration. However, the spatial patterns, the hemisphere differences in local summer monsoons and the associated mechanisms were not studied.

Here we focus on addressing the following questions: What are the spatial–temporal structures of the global scale local summer monsoon precipitation (GSMP) change on the timescales longer than a decade over the past millennium? Can we distinguish the forced response and internal variability of the GSMP? What are the driving mechanisms for the GSMP as well as the NHSMP and SHSMP? Are there any difference in the driving mechanisms between the NHSMP and SHSMP?

2 The model and validation

2.1 Model and experiments

The numerical model used for millennial integration is the ECHO-G coupled climate model (Legutke and Voss 1999), which consists of the spectral atmospheric model ECHAM4 (Roeckner et al. 1996) and the global ocean circulation model HOPE-G (Wolff et al. 1997). The description of the model physics and performance was given in details by Zorita et al. (2003). Here for convenience of the readers, we provide a brief summary. The model configuration has 19 vertical levels in the atmosphere and 20 levels in the ocean, and horizontal resolutions are approximately 3.75° (atmosphere) and 2.8° (ocean) in both latitudes and longitudes. The ocean model HOPE-G has a grid refinement in the tropical regions, where the meridional grid point separation reaches 0.5°. To enable the coupled model to sustain a simulated climate near to the real present day climate with minimal drift, both heat and fresh-water fluxes between atmosphere and ocean are modified by adding a constant (in time) field of adjustment with net-zero spatial average (Roeckner et al. 1996; Wolff et al. 1997).

Two millennial integrations with the ECHO-G model were conducted: a 1000-year control simulation (CTL) which was generated using fixed external (annually cycling) forcing set to the present-day values (Zorita et al. 2003) and a forced run, named ERIK, covering the period 1000–1990, which is forced by three external forcing factors (Gonzalez-Rouco et al. 2003; von Storch et al. 2004; Zorita et al. 2005): solar variability (Crowley 2000), greenhouse gas (GHG) concentrations in the atmosphere including CO2 and CH4 (Blunier et al. 1995; Etheridge et al. 1996) and the effective radiative effects from stratospheric volcanic aerosols (Crowley 2000) for the period 1000–1990 AD (Fig. 1). The volcanic forcing is parameterized in this simulation as a simple reduction of the annual mean solar constant, starting in the year with a volcanic eruption and usually lasting a couple of years, according to the reconstructions of volcanic aerosol forcing (Crowley 2000).

Fig. 1
figure 1

The external forcing used in ECHO-G forced run: the effective solar radiation forcing (insolation plus volcanic aerosol effect; solid line), the concentration of CH4 (long dashed line) and CO2 (dotted line)

The initial conditions of the ERIK simulation were taken from year 100 of the control run. Those initial conditions are, however, representative of present-day rather than pre-industrial climate and the experimental design therefore included a 30-year adjustment period during which the control run forcing was linearly reduced until it matched the forcing imposed around 1000 AD, followed by a 50-year period with fixed forcing to allow the model’s climate to readjust to the modified forcing. The ERIK simulation then proceeded from the conditions at 1000–1990 AD. Note that the uncertainty in the initial conditions, in turn, could potentially influence the relationship between applied forcing and simulated response in the first 100–200 years. This issue has been addressed by Zorita et al. (2007).

2.2 Validation of the model precipitation climatology

How the model responds to annual variation of the insolation is critical for model evaluation because we examine how the model responds to long-term variability in insolation forcing. The performance of ECHO-G in modeling the patterns and errors of the climatological annual mean and annual cycle of precipitation has been assessed by Liu et al. (2009). Their results indicate that the simulated precipitation climatology in ERIK run is comparable to those assimilated climatology in NCEP-2 reanalysis. The overall agreement adds confidence to our further analysis of the decadal-centennial precipitation variability using the outputs generated by the ERIK run. Despite these successes, the simulated annual mean precipitation has notable biases in East Asia monsoon region, subtropical South Pacific and South Atlantic convergence zones, and the Mexican-North American monsoon regions (Liu et al. 2009). These biases are also common in most of coupled models in simulating annual mean precipitation (Lee et al. 2010).

Figure 2 compares the monsoon domain defined by using CMAP (Climate prediction center Merged Analysis of Precipitation) climatology and by the climatology simulated in ERIK (forced) run. Evidently, the simulated monsoon precipitation domains are reasonably realistic. Major deficiency is the missing of the East Asian subtropical monsoon. Other minor differences are mainly due to different model (3.75 × 3.75°) and CMAP (2.5 × 2.5°) resolutions. Generally speaking, the model simulations on global scale are adequate for our study of the long-term modulations of the global monsoon under imposed natural and anthropogenic forcing. For convenience, in the following analysis we will use CMAP monsoon domain.

Fig. 2
figure 2

Comparison of the global monsoon precipitation domains obtained by using a CMAP (1979–2008) and b ERIK simulation output (1000–1990). The monsoon domain is the regions where the annual range (local summer mean minus winter mean) of precipitation rate exceeds 2.0 mm/day and the local summer precipitation exceeds 55 % of the annual total. Here local summer is defined as May through September (MJJAS) for the NH and November through next March (NDJFM) for the SH. The dotted regions represent dry regions where the local summer precipitation rate is below 1 mm/day

2.3 Statistical significance test for the difference of the regression coefficients

We used linear regression to estimate trends for given periods. Therefore test of the significance between two trends becomes a significance test of the differences between two regression line slopes, or regression coefficients. Denote the two slopes by b 1 and b 2. To test whether the difference between the two slopes is not due to sampling errors, we use the Student’s t test with

$$ t(n_{1} + n_{2} - 4) = \frac{{b_{1} - b_{2} }}{{S_{yx} \sqrt {\frac{1}{{S_{x1} }} + \frac{1}{{S_{x2} }}} }} , $$
(1)

Where S x1 and S x2 denote the sums of squares for x in the first and second group, respectively, and

$$ S_{yx} = \sqrt {\frac{{S_{y1} + S_{y2} - b_{1}^{2} S_{x1} - b_{2}^{2} S_{x2} }}{{n_{1} + n_{2} - 4}}} , $$
(2)

where S y1 and S y2 are the sums of squares for y in the first and second groups, and n 1 and n 2 are sample sizes for these two groups, respectively.

3 The forced and internal modes of the decadal GSMP

Can we distinguish forced and internal variability of the GSMP? The separation is possible if the forced response and internal variability have distinct spatial patterns and temporal behaviors. For this reason, we examine the principal modes of the GSMP variability by empirical orthogonal function (EOF) analysis of the decadal mean precipitation in monsoon domains over the past millennium. The first and second modes accounts for 9.5 and 5.4 % of the total variance, respectively, and they are significantly separated from each other and from the rest higher modes using the method of North et al. (1982).

Figure 3 shows the leading mode of GSMP variability and associated (regressed) SST and 850 hPa wind anomalies. The spatial structure is characterized by a nearly uniform increase (or decrease) of monsoon precipitation across all regional monsoon regions, which is in association with a global SST warming (or cooling). The global SST pattern features obviously warming (or cooling) in the mid-latitude oceans. But the warming (or cooling) signal is relative weak in the region from the southeastern to equatorial western-central Pacific. The corresponding principal component (PC1) shows a slow variation on centennial-millennial time scale with enhanced GSMP generally occurring in the MWP from 1000 to 1400 AD and the twentieth century (1900–1990 AD), along with a suppressed monsoon rainfall occurring in the LIA from 1450 to 1850 AD. Spectral analysis reveals a bicentennial peak significant at 90 % significance level (Fig. 3d). The sub-millennial variation and the bicentennial peak are primarily a response to the effective solar forcing because the SV forcing shows a similar sub-millennial variation and the spectrum of the SV forcing has a similar bicentennial peak (Liu et al. 2009). Of note is that the increase of GSMP during the twentieth century is unprecedented, resulting in a maximum that is significantly stronger than that during the MWP. Since the effective solar forcing in the MWP is about the same as in the twentieth century, the excessive warming in the twentieth century reflects the impact of the GHG forcing. Note also that this mode is totally absent in the CTL (forcing-free) run, indicating that it is externally forced by effective solar radiation (Solar-Volcanic, SV) and GHG forcing.

Fig. 3
figure 3

The leading EOF mode of the GSMP variation over the past millennium: a spatial structure, b the corresponding principal component (PC); the smoothed curve is 3-decade running mean, c regressed SST and 850 hPa wind anomalies, shown are only wind vectors that are significant above 90 % confidence level; and d the corresponding PC spectra. The fraction on the top of panel a denotes the fractional variance. The data were decadal mean time series derived from ECHO-G forced (ERIK) simulation. The correlation matrix was used for EOF analysis

The second mode does not show a cohesive pattern of variability (Fig. 4). The NH monsoon rainfall tends to be in phase while the SH monsoon rainfall tends to be out of phase between the eastern and western hemisphere. The corresponding principal component (PC2) has triple spectral peaks on 2–7 decades (Fig. 4d). The increasing precipitation pattern is associated with a global SST rising with a maximum warming in the equatorial western-central Pacific (5°S-5°N, 150°E-120°W) and a weak cooling in the equatorial eastern Pacific as well as in the northern and southern oceans. This concentration of SST variability in the equatorial central Pacific is quite different from the SST pattern associated with the EOF1 where midlatitude warms more than the equatorial regions. What is the origin of the EOF 2 mode?

Fig. 4
figure 4

The same as in Fig. 3 except for the second EOF mode of the decadal mean GSMP

To address this question, we examine the control (CTL) run in which the external forcing is fixed. Figure 5 shows the leading EOF mode derived from the control experiment. Because the forcing is fixed in the CTL run, the leading mode must arise from internal feedback processes within the coupled climate system. We note that the spatial structure of the EOF 2 mode in the forced run (Fig. 4a) is very similar to that of the leading mode in the free run (Fig. 5a). The regressed SST patterns are also very similar (figure not shown). In general, the power spectra are also similar. This suggests that the EOF 2 mode in the forced run is essentially an internal variability mode. The high degree of similarity in the internal modes between the forced and free runs suggests that the internal variability mode is not significantly changed by the presence of external forcing.

Fig. 5
figure 5

The leading EOF mode of the GSMP variation over the past millennium derived from the ECHO-G CTL (forcing-free) run: a spatial structure, b the corresponding principal component with 3-decade running means, and c the corresponding spectrum. The number on the top of panel a denotes the fractional variance

4 Differences between the NHSMP and SHSMP variations

In this section, we investigate the changes of the total amount of monsoon precipitation in the NH, SH, and the globe, i.e., NHSMP, SHSMP and GSMP in response to the external forcing, which drives both the global mean 2m-air temperature (GMT) and GSMP (Fig. 6b).

Fig. 6
figure 6

a The four forcing (decadal mean) time series: solar (red, W/m2), volcanic (blue, W/m2), effective solar forcing (black, W/m2), and CO2 concentration (purple, ppm); (b) The global mean 2m-air temperature (GMT, red, °C), GSMP (black in (i), mm/day), and NHSMP (black in (ii), mm/day), and SHSMP (black in (iii), mm/day), respectively. The smoothed curves are 3-decade running means. The correlation coefficients with GMT are shown at the lower-right corner of each panel (the correlation coefficients are calculated using 3-decade running mean time series)

Figure 6b shows that the change of GSMP over the past millennium follows that of the GMT closely (r = 0.87). This in-phase variation occurs primarily on the centennial-millennial time scale because the interdecadal variations are dominated by the internal feedback process which does not contribute significantly to global mean quantities. The GMT decreases from 1140 to 1700 AD and then increases from 1700 to 1990 AD. Similarly, the GSMP decreases also roughly from 1140 to 1700 AD and rises from 1700 to 1990 AD. The rapid increase of GMT since around 1850 seems to be associated with the rapid increase of CO2 concentration.

Note that the global monsoon is the dominant mode of the annual variation of the tropical precipitation and circulation (Wang and Ding 2008) and the GSMP intensity measures change of the annual cycle of the coupled climate system. Therefore, the in-phase relationship between the GMT and GSMP indicates that global warming is accompanied by amplification of the annual cycle of the Earth’s climate system. Furthermore, the GMT measures the annual mean climate but the GSMP measures the annual cycle of the climate system. Thus, the GSMP provides a sensible and complementary measure of global climate change in terms of a key hydrological variable.

It is of interest to examine the rate of change in NHSMP and SHSMP as functions of GMT (precipitation change rate per K). We found that the ratios of NHSMP/GMT and SHSMP/GMT tend to differ in response to the external SV forcing. As shown in Fig. 7, during the preindustrial period (1000–1850 AD) during which the SV forcing dominates and GHG forcing has negligible change, the ratio NHSMP/GMT is about 2.5 % K−1, while the ratio SHSMP/GMT is about 4.1 %K−1. The difference is statistically significant at 95 % confidence level by the Student’s t test. The result means that the SHSMP is more sensitive to the SV forcing than the NHSMP.

Fig. 7
figure 7

Scatter diagram showing NHSMP and SHSMP rates as functions of global mean 2m-air temperature for the preindustrial period 1000–1850. The difference between the two slopes is significant at 95 % confidence level, suggesting that the SHSMP is more sensitive to the SV forcing than the NHSMP

On the other hand, we found that the NHSMP responds to GHG forcing more sensitively than the SHSMP. Figure 8 shows the NHSMP and SHSMP rates as functions of GMT for the preindustrial period (1000–1850 AD) and the industrial period (1850–1990 AD) during which the GHG forcing increases sharply. During the preindustrial period, the ratio NHSMP/GMT is about 2.5 % K−1, while it is about 4.3 %K−1 in the industrial period (Fig. 8a). The difference between these two is statistically significant at 95 % confidence level. Since the difference is due to the presence of the GHG forcing in the industrial period, we infer that the NHSMP responds to the GHG forcing more sensitively than to the SV forcing. On the other hand, the ratios SHSMP/GMT in the preindustrial and industrial periods are about 4.1 and 4.5 %K−1, respectively (Fig. 8b), and the difference between these two is statistically insignificant, suggesting that the SHSMP is insensitive to the GHG forcing. So the NHSMP responds to GHG forcing more sensitively than the SHSMP does.

Fig. 8
figure 8

Scatter diagram showing a NHSMP and b SHSMP rates as functions of global mean 2m-air temperature for the preindustrial period 1000–1850 (blue squares) and the industrial period 1850–1990 (red circles). The difference between the two slopes is significant at 95 % confidence level for NHSMP in a, but not significant for SHSMP in b

Why is the NHSMP more sensitive to GHG than to SV forcing? It is conceivable that the GHG forcing creates larger land–ocean thermal contrast than SV forcing. This has been tested by calculation of the thermal contrast changes in preindustrial period (MWP minus LIA) and the industrial period (Present minus LIA). The results show that the NH land–ocean thermal contrast is 0.18 °C for MWP minus LIA and 0.48 °C for Present minus LIA, indicating that GHG forcing is indeed more effective than SV forcing in creating land–ocean thermal contrast, and thus favors increased NHSMP.

5 Factors controlling the NHSMP, SHSMP, and GSMP

What controls the centennial-millennial variation of the GSMP? Are there any differences between the NHSMP and SHSMP? Figures 9, 10, and 11 provide answers to these questions by presenting the linkages of the NHSMP, SHSMP and GSMP with the lower boundary conditions, such as the sea level pressure (SLP) and 2m-air temperature (T2m), both reflecting the land surface thermal conditions and SST conditions faithfully.

Fig. 9
figure 9

Factors controlling NHSMP. The decadal mean time series of a NHSMP (mm/day), b land–ocean thermal contrast: T2m (land) minus T2m (ocean) between the equator to 60°N in MJJAS (°C), and c hemispherical thermal contrast: T2m (NH, 0–60°N) minus T2m (SH, 0–40°S) in MJJAS (°C). The thick lines denote 3-decade running means. The correlation coefficients with NHSMP using 3-decade running mean time series are shown at the lower-right corner of each panel

Fig. 10
figure 10

Factors controlling SHSMP. The decadal mean time series of a SHSMP (mm/day), b the east–west mass contrast: SLP in southeastern Pacific (0–40°S, 180–70°W) minus SLP in Indian Ocean (20°S–20°N, 40°E–120°E) during NDJFM (hPa), and c the SH subtropical high strength: SLP averaged between 20°S and 40°S (hPa). The thick lines denote 3-decade running means. The correlation coefficients with SHSMP using 3-decade running mean time series are shown at the lower-right corner of each panel

Fig. 11
figure 11

Factors controlling GSMP. The decadal mean time series of a GSMP, i.e., NHSMP + SHSMP (mm/day), b east–west mass contrast: annual mean SLP in southeastern Pacific (0–40°S, 180–70°W) minus SLP in Indian Ocean (20°S–20°N, 40°E–120°E) (hPa), c global land–ocean thermal contrast: annual mean T2m (land) minus T2m (ocean) between 40°S and 60°N in °C, d the SH subtropical high strength: annual mean SLP (40°S–20°S). The thick lines denote 3-decade running means. The correlation coefficients with GSMP using 3-decade running mean time series are shown at the lower-right corner of each panel

Figure 9 indicates that a stronger NHSMP is related to (a) an increased land–ocean thermal contrast (warm land and cold ocean in the NH), and (b) an increased hemispheric thermal contrast (warm NH and cold SH). Increased land–ocean thermal contrast can enhance monsoon lows and the associated moisture convergence (Liu et al. 2009); and warmer NH generates cross-equatorial pressure gradients that drive low-level cross-equatorial flows from SH to NH, again strengthening the NH monsoon (Loschnigg and Webster 2000). The combined effects of the land–ocean and NH-SH thermal contrasts yield a strong control of the NHSMP variability on decadal-centennial time scales.

Different from the NHSMP, the SHSMP is not significantly related to the land–ocean and hemispheric thermal contrasts (figure not shown). Instead, it is associated with (a) an east–west thermal contrast-induced SLP difference between the southeastern Pacific and tropical Indian Ocean, and (b) the circumglobal SH subtropical high strength measured by the SLP averaged between 20°S and 40°S (the ridge lines of the SH subtropical high locate at about 30°S) (Fig. 10). The SHSMP is enhanced when pressure rises in the southeastern Pacific and drops in the tropical Indian Ocean. The increased east–west SLP gradient induces a westward air flow and moisture transport, which enhances moisture convergence over the southern African-southwest Indian Ocean and the Australian monsoon region. The enhanced SH subtropical high increases the trade winds in SH oceans, also strengthening the convergence of air mass and moisture in the SH monsoon low pressure regions.

The GSMP is the sum of the NHSMP and SHSMP, so it is affected by the common factors that control the NH and SH summer monsoons. Analysis shows that the GSMP is most closely associated with (a) the east–west thermal contrast-induced SLP difference between the southeastern Pacific and tropical Indian Ocean (r = 0.88), (b) the land–ocean thermal contrast over the globe between 40°S and 60°N (r = 0.84), and (c) the circumglobal SH subtropical high strength between 40°S and 20°S (r = 0.88) (Fig. 11). The enhanced “southeastern Pacific cooling-Indian Ocean warming” induces rising SLP in the eastern Pacific, which enhances the trade winds, transporting and converging moisture into the eastern hemisphere monsoon regions, which not only favor enhancement of the SH monsoon (Fig. 10) but also NH monsoon through enhancing cross equatorial flows. The increased land–ocean thermal contrast enhances monsoon low and associated moisture convergence, and primarily favors the NH summer monsoon precipitation. The enhanced SH subtropical high strengthens not only the SHSMP (Fig. 10) but also the NHSMP by increasing the pressure gradient between the NH and SH. The combined effects of the east–west and the land–ocean thermal contrasts, as well as the SH subtropical high strength result in a strong control of the GSMP.

6 Concluding remarks

The GSMP measures the annual range of precipitation averaged in the GM domain, thereby providing a useful parameter to quantify the change in the annual cycle of the earth climate. This new measure also provides important information about the change in the global monsoon precipitation, which is a climate variable that is far more relevant for food production and water supply than the mean temperature change. Note that the GMT measures change of annual mean climate in terms of temperature, whereas GSMP measures change of the annual cycle in terms of precipitation. They are complementary in nature. The linkage between the GMT and GSMP indicates the close relationship between the global warming and amplification of the annual cycle of the climate system.

We have shown that forced response of the GSMP has a distinct spatial–temporal structure from that due to internal feedback processes. Two prominent patterns of GSMP variability are identified. The leading pattern has a dominant centennial-millennial variation and features a nearly uniform increase of monsoon precipitation across all regional monsoons. This pattern is a forced response to the changes in effective solar-volcanic (SV) radiation and GHG concentration. The second pattern is associated with a multi-decadal oscillation in the central Pacific sea surface temperature (SST), representing an internal feedback mode.

An important finding is that the NHSMP and SHSMP have notable differences in terms of their driving mechanisms. First, the NHSMP responds to GHG forcing more sensitively than the SV forcing, while the SHSMP responds to the natural solar-volcanic radiative forcing more sensitively than the NHSMP does. Second, the enhanced NHSMP is primarily driven by warm land-cold ocean in NH and warm NH-cold SH, while the enhanced SHSMP is primarily caused by enhanced east–west thermal contrast over the SH Indo-Pacific warm pool and the SH subtropical high strength. The GSMP is driven by the factors that commonly control both the NH and SH summer monsoons, including the east–west thermal contrast between the Southeast Pacific and tropical Indian Ocean, the global land–ocean thermal contrast, and the circumglobal SH subtropical high strength.

These results obtained from the present pilot study carried out by using ECHO-G millennial runs should be compared with available proxy data and tested by using other models and by analyzing model outputs on other time scales. These include examination of PMIP-3 ensemble model outputs for past millennium, last glacial maximum and mid-Holocene periods.

The physical mechanisms leant from this study will also add understanding of future change of the global monsoon. For instance, the climate models’ future projections under increasing GHG forcing display an increased land–ocean thermal contrast and the contrast between warmer NH and colder SH. These two factors in principle may favor enhancement of the NHSMP, and the GSMP as well.