Introduction

High-quality social development with effective disaster prevention, logistics-planning, agriculture production, and decision-making requires accurate sub-seasonal rainfall prediction over land monsoon regions where about two-thirds of the world’s population live1. Tropical intraseasonal variabilities (ISVs) and the associated teleconnections are dominant sources of sub-seasonal predictability2,3,4.

Madden-Julian Oscillation (MJO)5,6 is a major mode of tropical ISV2,7. During boreal summer, MJO’s eastward propagation significantly weakened, meanwhile pronounced northward and northeastward propagation prevails in the Northern Hemispheric (NH) monsoon regions8,9,10,11,12,13. The MJO’s impacts on regional monsoon rainfall have been extensively investigated. These include South Asian (SA)14,15,16, East Asian (EA)17,18,19, North American (NAM)20,21,22, northern African (NAF)23,24,25, Australian (AU)26,27,28, South American (SAM)29,30,31,32,33,34, and South African (SAF)35,36,37 monsoons.

Prediction of weekly-mean land monsoon rainfall at a lead of 2–8 weeks has been a major effort yet particularly challenging. The current evaluation of sub-seasonal prediction skill and predictability estimates has primarily focused on the leading modes of the MJO or boreal summer intraseasonal oscillation3,38,39,40,41. The state-of-the-art European Centre for Medium-Range Weather Forecast (ECMWF) model can predict the leading modes of the MJO up to 40 days in advance in terms of the Real-time Multivariate MJO (RMM) index42,43. However, for up to four weeks, the MJO’s rainfall prediction skill is confined to a narrow equatorial belt, and the rainfall prediction skills drop rapidly beyond two weeks in vast off-equatorial monsoon regions44,45. User community concerns with weekly rainfall prediction rather than the MJO index. The reason behind the limited sub-seasonal predictability of monsoon rainfall is the forefront challenge faced by the sub-seasonal prediction community.

This further raises the question of why the weekly-mean monsoon rainfall prediction skill declines much faster than the MJO index skill in the prediction models. The MJO index measures the equatorial region’s low-frequency (LF, ~20–70 day) ISV. In most off-equatorial monsoon regions, the MJO index only accounts for a fraction of the total intraseasonal rainfall variability46. We suspect that the Week 2 to Week 3 monsoon rainfall is mainly dominated by the high-frequency (HF, ~8–20 day) ISV, not the MJO. We might not have fully recognized the importance of HFISV in sub-seasonal prediction and predictability. The HFISVs may originate from convectively coupled equatorial Rossby waves, named quasi-biweekly oscillation47, and from mid-latitude wave trains due to eddy-mean flow interaction, eddy–eddy interaction, and others48,49,50. The HFISVs in various monsoon regions have been documented. Over SA, HFISV was dominated by the westward-propagating quasi-biweekly oscillation16,47,51,52,53. The significant HFISV of EA precipitation was linked to preceding mid-latitude wave trains and tropical quasi-biweekly oscillation18,53,54,55,56. HFISVs were also reported over NAM57, NAF23,24,58, AU27, and SAM29,33,59,60,61. However, the predictability sources of HFISV and the processes by which the equatorial disturbances and mid-latitude processes affect the HFISV of different monsoons are not well understood.

The ISVs have been investigated primarily in individual monsoon regions, as mentioned above, because each regional monsoon has an indigenous land-ocean configuration and involves different atmosphere–ocean–land interaction processes and convective and mesoscale systems62. On the other hand, an across-the-board synthesis of the common and distinct features among the regional monsoons’ ISVs is persuasive for an in-depth understanding, simulation, and prediction of the ISVs. The eastward-propagating planetary-scale MJO circulation system and mid-latitude disturbances may significantly regulate and coordinate the regional ISVs by changing monsoon circulations. These coordination and regional influences entail a global perspective. An overview of the ISVs of the land monsoon rainfall from a global perspective is advantageous for a deeper understanding of MJO–monsoon interaction, different characteristics and sources of regional ISV, and especially for improving the sub-seasonal prediction of regional land monsoon rainfall.

Our particular interests also include examining recent changes in regional monsoon ISV intensity under global warming. The seasonal-mean rainfall in the NH monsoon regions has increased considerably since 197963, and the MJO residence time was observed to increase over the Indo-Pacific Maritime Continent by 5–6 days during the last four decades due to the warm pool expansion64. These changes of mean state and MJO might have influenced the land monsoon ISV. The regional ISV was found to be enhanced since 1993 over South China65. The changes in both HF and LFISV intensities in most land monsoon regions, however, remain unknown.

This study reviews and investigates the rainfall ISVs among individual land monsoons, including different leading modes of variability, their origins, tracks of HFISV disturbances, and recent changes in ISV intensity. We also try to explain the divergent sub-seasonal prediction skills over different land monsoon regions in current sub-seasonal to seasonal (S2S) prediction models. The results are expected to deepen understanding of the dynamics of the monsoon ISV and improve ISV’s simulation and sub-seasonal prediction of land monsoon rainfall.

Results

Relative importance of HF and LFISVs in summer monsoon regions

The maximum intraseasonal rainfall anomalies tend to anchor over the global land monsoon regions (Fig. 1a). The average ISV intensity of the global land monsoon rainfall, measured by the standard deviation of the 8–70-day filtered precipitation anomaly, is 4.6 mm day−1, ranging from 3.6 mm day−1 (NAF monsoon) to 5.9 mm day−1 (AU monsoon) (Table 1). The ISV intensity is highly correlated with the corresponding seasonal-mean monsoon precipitation (Supplementary Fig. 1). The ISVs contribute to one-quarter to one-third of the daily precipitation variance over various land monsoon regions (Fig. 1b), ranging from 20.5% (NAF monsoon) to 32.4% (AU monsoon) (Table 1). The ISVs have relatively large contributions over AU, South and East Asian, Mexican, and southeastern SAM monsoons but low contributions over NAF, northwestern SAF, and western SAM monsoons.

Fig. 1: Climatological ISV activity in monsoon seasons during 1979–2018.
figure 1

a is summer (MJJAS for NH and NDJFM for SH) ISV intensity (mm day−1) defined by the standard deviation of 8–70-day bandpass-filtered precipitation anomalies. The global land monsoon regions are outlined by red lines, and the East Asian (EA) and South Asian (SA) regions are separated by the gray line. A monsoon region is defined as the region where the local summer-minus-winter precipitation exceeds 300 mm, and the summer precipitation exceeds 55% of the annual total103. b shows the 8–70-day ISV variance contribution to the total daily precipitation variance. The contribution is calculated by the ratio of the averaged power spectrum in the intraseasonal band to the total power. The area-conserving format is used in which the variance is proportional to the area of the logarithm of frequency versus power times frequency.

Table 1 Averaged ISV intensity (measured by the standard deviation of 8–70-day bandpass-filtered precipitation anomalies shown in Fig. 1) and fractional variance contribution of ISV to the total daily variance for each regional land summer monsoon.

The regionally averaged land monsoon precipitation spectra show a significant 8–15-day HFISV peak in NH monsoons and a significant 10–20-day HFISV in SH monsoons (Fig. 2). The precipitation spectra also show significant LF peaks in AU (20–70 days) and in SA and SAM (20–45 days). No significant signal is observed beyond 70 days. Thus, we define the HF component with a time scale of 8–20 days (~2–3 weeks) and the LF component, 20–70 days (~4–9 weeks).

Fig. 2: Power spectra of ISVs of individual summer monsoons.
figure 2

Power spectra of the 3-day running mean precipitation anomalies (mm day−1) averaged over each monsoon region for a NH and b SH during the summers of 1979–2018. Daily climatology is removed to obtain the anomalies. Each monsoon precipitation series is normalized before the spectrum analysis. Dashed lines denote the corresponding 90% confidence level.

The relative contributions of these HF and LFISVs to the total ISV variance are defined by the ratios of their associated power spectra to the total (see Methods), based on the area-conserving format of power spectra (Supplementary Fig. 2). The 8–20-day ISV accounts for 53–70% of the total 8–70-day ISV variance in each regional land monsoon, whereas the 20–70-day ISV accounts for 30–47% (Table 2). The EA monsoon exhibits the largest 8–20-day ISV contribution (70.2%), and the AU monsoon has the highest 20–70-day ISV contribution (47.3%). Over AU, SA, and SAM, the HFISV contribution (53–60%) is slightly higher than LFISV (40–47%), but over EA, NAF, NAM, and SAF, the HFISV accounts for about 66–70% of the total ISV variance, dominating the ISV. The results suggest that MJO has a more significant influence on AU, SA, and SAM monsoons but less on the other four monsoon regions. The proximity of the AU monsoon to the MJO convective activity during austral summer explains why the AU monsoon has a strong LFISV component. During boreal summer, the equatorial Indian Ocean remains an MJO activity center. The prominent monsoon vertical wind shear and the northward shift of the maximum moist static energy zone favor the equatorial convective anomaly moving northward to affect SA ISV. Over the SAM, the MJO-induced mid-latitude teleconnection likely has a large amplitude34,66,67, affecting SAM through the northward propagating cyclonic circulation anomalies. Due to the lower contribution of MJO to monsoon ISVs in EA, NAM, NAF, and SAF, it is interesting to find out whether the rainfall sub-seasonal predictability in these monsoon regions is also significantly lower than that in AU, SA, and SAM.

Table 2 Fractional variance of the 8–20-day and 20–70-day ISVs against the total 8–70-day ISV for averaged summer monsoon precipitation over each monsoon region.

The Maritime Continent is part of the monsoon region in terms of its annual reversal of the prevailing wind68,69. The seasonal-mean precipitation and 8–70-day ISV intensity averaged over the Maritime Continent (10oS–5oN; 100o–150oE) are 12.5 mm day−1 and 5.8 mm day−1, respectively, during its wet season from November to March (Supplementary Fig. 3a). The Maritime Continent rainfall only shows a significant 8–20-day HFISV component during the wet season (Supplementary Fig. 3b), which contributes to 58.3% of the total 8–70-day variance.

Origin and propagation of 8–20-day ISV in each monsoon region

We first identify the leading spatial pattern and activity center of the HFISV of summer precipitation in each monsoon region. Figure 3a shows that the dominant EOF modes of the HFISVs of land monsoon precipitation exhibit a uniform structure in the SH monsoon and NAF monsoon regions but a north-south dipolar structure over the SA, EA, and NAM monsoon regions. The SA monsoon is characterized by wet northern India-Bangladesh and dry Central India. The EA monsoon shows a contrast between a wet Southeast China and a dry Indo-China Peninsula. The NAM monsoon features a wet Venezuela and a dry Mexico. Conceivably, the North-South dipolar structures occur in the monsoon regions with a larger meridional extent, including both tropics and subtropics, so that one center is in the tropical monsoon trough and the other in the subtropical convergence zone (poleward of 20o latitude). The uniform precipitation anomalies in the NAF, AU, SAM, and SAF monsoons are centered in Nigeria, northern Australia, Brazil, and Mozambique-Madagascar, respectively. These leading modes can explain about 9% of the 8–20-day precipitation variance, ranging from 11.4% in EA to 7.2% in SAF. Note that the second EOF modes have a comparable fractional variance to the first ones (Supplementary Table 1). Together, the two leading modes explain about 15% of the variance.

Fig. 3: Dominant ISV modes in individual summer monsoon regions.
figure 3

The spatial patterns of the first EOF modes of summer sub-seasonal precipitation anomalies (shading; mm day−1) for a 8–20-day ISV and b 20–70-day ISV over each monsoon region during the summers of 1979–2018. The EA and SA regions are separated by the gray line.

To investigate where the peak phases of the leading 8–20-day ISV modes originate, we use the sequential leading maps of regressed OLR anomalies to trace their origins and propagation pathways from the tropics, and regressed wind anomalies to trace the leading signals from both the tropics and mid-latitudes (Figs. 4, 5), following many previous works47,55,70. On Day 0, the regressed OLR anomalies are consistent with the precipitation anomalies of the leading 8–20-day ISV modes shown in Fig. 3a.

Fig. 4: Origins and propagations of 8–20-day ISVs in NH summer monsoon regions.
figure 4

Shown are the lead-lag regression maps of 8–20-day bandpass-filtered OLR anomalies (shading; W m−2) and wind anomalies (vector; m s−1) onto the PC1s during the summers of 1979–2018 for a SA, b EA, c North American (NAM), and d northern African (NAF) monsoons. The wind anomalies shown in a, b, and c are at 850 hPa, while those in d are at 200 hPa. Only OLR and wind anomalies significant at the 90% confidence level are shown. Letter A tracks the center of the anomalous anticyclone associated with the monsoon ISV, C the cyclone, S the southerly, N the northerly, E the easterly, and W the westerly.

Fig. 5: Origins of 8–20-day ISVs of SH summer monsoons.
figure 5

Same as Fig. 4, except for 850-hPa wind anomalies in SH: a Australian (AU), b South American (SAM), and c South African (SAF) monsoons. The red line denotes the trough center at 850 hPa.

Over SA (Fig. 4a), the low-level anticyclonic anomaly that couples with the suppressed convection emerges over the Philippine Sea around 135oE on Day -8. It then moves westward along 10–20oN, finally forming a suppressed rain band stretching from Central India to the Bay of Bengal on Day 0. The southwesterly anomalies to the north of the anticyclonic circulation strengthen moisture transport and convergence, leading to enhanced convection over northern India-Bangladesh. The zonally elongated anticyclone stretches ~50o in longitude, and the westward phase speed is about 7.2 m s−1. This feature is consistent with the previous findings16,51,52. An apparent upper-tropospheric wave train also moves westward along 40oN from eastern China to the Iranian Plateau against the background westerlies (Supplementary Fig. 4a), suggesting that the upper-level circulation anomalies are likely a response to the convective heating rather than an extratropical forcing. This strong upper-tropospheric wave train indicates that the SA precipitation can significantly modulate the mid-latitude circumglobal teleconnection71, “silk road” teleconnection72,73, and the Tibetan Plateau74 on the HF intraseasonal time scale.

In EA (Fig. 4b), an anticyclonic anomaly emerges from the western equatorial Pacific around 140oE on Day -8; it propagates northwestward and is accompanied by suppressed convection, reaching the South China Sea and Indo-China Peninsula on Day 0. The westward phase speed is about 3.5 m s−1, consistent with previous result75. This northwestward propagation pathway is similar to the 8–10-day disturbances documented by Lau and Lau76. In theory, moist Rossby waves can emanate from the equatorial convective anomalies when the equatorial intraseasonal convective anomalies decay, passing through the western Pacific10. They then move northwestward under the influence of background monsoon easterly vertical wind shear shown in Fig. 4 in Wang and Xie10. Note that from Day -2 to Day 0 convection over Southeast China is enhanced by the convergence of southwesterly wind anomalies to the northwest of this tropical anticyclonic anomaly, forming a dipolar structure. In addition to the tropical origin mentioned above, a lower-tropospheric cyclonic wind anomaly emerges over Lake Baikal at Day -8 and propagates southeastward (Fig. 4b). The southwest part of this cyclonic anomaly triggers convection over eastern China on Day -2, and the convection peaks on Day 0. Meanwhile, the major body of the cyclonic anomaly propagates into the North Pacific. The southeastward propagation is also obvious in the upper troposphere (Supplementary Fig. 4b). The evidence here seems to confirm previous findings concerning the influence of the southeastward propagation of mid-latitude wave trains on the quasi-biweekly variability of EA summer monsoon54,55,56,77.

Over the NAM and NAF monsoon regions, the OLR anomalies are relatively weak (Fig. 4c, d). However, their precipitation anomalies are not (Fig. 3a), suggesting that shallower convection may prevail over these two regions, a feature similar to the corresponding seasonal-mean precipitation78,79,80. For the NAM monsoon, an incipient low-level cross-equatorial southerly anomaly emerges in the central equatorial Atlantic on Day -6, which couples with enhanced convection and propagates westward, reaching a peak over Venezuela on Day 0 (Fig. 4c). The coupling of the enhanced convection with an upper-level cross-equatorial northerly anomaly is more obvious (Supplementary Fig. 4c). Meanwhile, a suppressed convection propagates from Venezuela on Day -8 to the Gulf of Mexico on Day 0, initially coupling with a cross-equatorial northerly and later with an anticyclonic anomaly. Along the equatorial Atlantic, the convective anomaly coupled with the cross-equatorial flow resembles a mixed Rossby-gravity wave81. We suggest that the 8–20-day ISV in the NAM region originates from an equatorial mixed Rossby-gravity wave and gradually transforms into a convectively coupled Rossby wave. No significant mid-latitude wave train is observed (Supplementary Fig. 4c).

The HFISV over NAF has a weak convective anomaly, which appears to be excited by upper-level divergence associated with a fast eastward-propagating equatorial Kelvin wave (Fig. 4d). The upper-level Kelvin wave easterly anomalies emerge over the eastern equatorial Pacific around 120oW on Day -6, propagate eastward, and reach the eastern Atlantic around 10oW on Day 0, with a phase speed of about 24 m s−1. The convection over NAF is then enhanced by the divergent flows associated with the easterly anomalies. The low-tropospheric westerly anomalies related to this Kelvin wave are relatively weak before reaching the Atlantic (Supplementary Fig. 4d).

Unlike in the NH, the HFISVs in all three SH monsoons are associated with significant preceding mid-latitude wave trains in the lower troposphere (Fig. 5) and upper troposphere (Supplementary Fig. 5). Over AU, the burst of monsoon convection on Day -2 is triggered by a low-pressure trough associated with a cyclonic anomaly over Central Australia. This cyclonic anomaly originates from the southeastern Indian Ocean on Day -8 and propagates eastward to Southeast Australia on Day 0 (Fig. 5a). In SAM, a cyclonic anomaly emerges to the west of Argentina along 40oS on Day -8, which propagates northeastward, triggering the convection over SAM on Day -6 (Fig. 5b). The convection anomaly couples with the cyclonic anomaly and matures on Day 0. The northeastward propagations are also associated with the sub-monthly oscillation of rainfall associated with the South Atlantic Convergence Zone59,60,82. The westward propagation to the east of SAM along 30oS seems to be a response to the HFISV of SAM. The role of mid-latitude wave trains on the ISV is also observed for the SAF monsoon (Fig. 5c). A mid-latitude low-level cyclonic anomaly propagates from the area over the Southeast Atlantic on Day -8 to the region over the southwestern Indian Ocean on Day 0. The cyclonic anomaly triggers the convection in Zimbabwe on Day -4 and then propagates northeastward and matures near Madagascar on Day 0. The corresponding propagations of the upper-tropospheric cyclonic anomalies for the three SH monsoons illustrate the mid-latitude wave trains’ equivalent barotropic structures (Fig. 5 and Supplementary Fig. 5). However, once coupled with convection in the monsoon region, the barotropic structure transforms into a baroclinic structure.

In summary, the ISVs in Asian monsoons are related to the preceding tropical quasi-biweekly oscillation originating from the western equatorial Pacific. The EA monsoon is additionally associated with southeastward-propagating mid-latitude wave trains. The HFISVs over NAM and NAF are triggered by the convectively coupled mixed Rossby-gravity and Kelvin waves, respectively. Besides, the westward propagation of HFISV from Venezuela to Mexico is accompanied by a transition from a mixed Rossby-gravity wave to an equatorial Rossby wave. All three SH monsoons are related to preceding mid-latitude wave trains in the SH.

The different origins of HFISVs are likely determined by monsoon locations and background circulations, as discussed over EA by Liu, et al.55. The SH monsoon HFISVs are triggered by the mid-latitude wave trains because of their proximity to the westerly jet stream (Supplementary Fig. 6). Likewise, the EA monsoon is also exposed to the mid-latitude wave train’s impact. On the other hand, NAM and NAF are close to the equator. Thus, their HFISVs originate from the convectively coupled equatorial waves.

MJO coordination of 20–70-day ISV in each monsoon region

The leading modes of LF (20–70-day) ISVs have similar spatial patterns as the corresponding HF (8–20-day) ISVs for each monsoon (Fig. 3). Nevertheless, minor differences exist. In SA, the Central Indian precipitation anomaly is stronger than that in northern India-Bangladesh (Fig. 3b). In NAM, the Mexican precipitation anomaly is stronger than the Venezuela anomaly. The similarity between the spatial patterns of the HF and LFISVs implies that the mean regional monsoon circulation might have similar controls on the HF and LFISVs.

In SA, AU, and SAM regions, where the LFISVs have significant contributions to the total ISV (Table 2), the first EOF mode explains 22.6, 18.5, and 17.0% of the total 20–70-day variance, more than double the fractional variance of the corresponding second modes (Supplementary Table 1). It suggests that the monsoon responses to the MJO over SA, AU, and SAM are geographically-locked amplification. This type of geographically-locked amplification has also been observed before for SA83 and for SAM84,85, which means that the mean regional monsoon circulation has significant control on the LFISVs. Conversely, in other monsoon regions where the HF dominates the ISV, the first two modes of the LFISVs have comparable fractional variances, suggesting a propagating mode.

The 20–70-day ISV in Central India is related to the northward propagation of the boreal summer ISO from the equatorial Indian Ocean (Supplementary Fig. 7a). The Southeast China precipitation is enhanced by the southwesterly wind anomaly of the western North Pacific suppressed ISO that originates from the Indian Ocean (Supplementary Fig. 7b), consistent with the previous works8,52,86. In SH, the AU monsoon 20–70-day ISV is dominantly affected by the MJO convection over the Arafura Sea (Supplementary Fig. 7c), as demonstrated by previous works28,87. The MJO modulates other monsoon regions primarily by exciting Kelvin waves. Significant preceding MJO convection anomalies are observed over the Warm Pool region. The MJO-excited Kelvin waves, represented by upper-level easterly anomalies, propagate eastward along the equator and enhance the convection over the Gulf of Mexico, Nigeria, Brazil, and Zambia-Madagascar (Supplementary Fig. 8).

Significant preceding mid-latitude wave trains are observed in SAM. The wave train associated with SAM originates from the ocean to the east of Australia and is tied to the MJO convection anomaly over the tropical western Pacific (Supplementary Fig. 8c). This Indo-Pacific MJO’s impact on SAM through the mid-latitude teleconnection has also been observed previously34,37,59,60,67,88,89,90. Since the group speed of the teleconnection is quite fast, and the MJO’s signal can reach America in about 1 week91,92 so that the SAM only lags the Indo-Pacific MJO by 10 days. It would be interesting to explore how much these mid-latitude LF wave trains are independent of the tropical MJO forcing in the future, as discussed by previous works93,94.

In summary, the 20–70-day ISVs of the Asian-Australian monsoon are directly affected by the MJO propagation. The Indo-Pacific MJO influences the American and African monsoons through its eastward-propagating Kelvin component and also has robust impacts on the SAM through its mid-latitude teleconnection in the SH.

Figure 6 summarizes the relationship between the MJO peak wet phase over the central equatorial Indian Ocean (or dry phase over the tropical western Pacific) and each region’s leading mode of 20–70-day ISV. All these correlations are significant, except for African monsoons. MJO peak phase has the strongest linkages with AU, SA, and SAM. The MJO peak wet phase over the Indian Ocean leads the AU wet phase by 14 days with the highest r = 0.3 (p < 0.01) and the SA wet phase by 12 days with r = 0.26 (p < 0.01), respectively. The MJO dry phase over the western equatorial Pacific, which concurs with the equatorial Indian Ocean wet phase, leads the SAM dry phase by ten days with r = 0.25 (p < 0.01). The high correlations in these regions are consistent with their higher fractional variance of the ISV to daily variance (Fig. 1b) and higher contribution of the LFISV to the total ISV than in other regions (Table 2). In addition, the enhanced MJO convection over the central equatorial Indian Ocean tends to lead wet Southeast Asia and dry South China by 22 days. The dry anomaly over the tropical western Pacific leads the dry phase by 15 days over Mexico and by ten days over Brazil. These lead-lag correlations can be potentially used for the sub-seasonal prediction of each regional monsoon.

Fig. 6: Role of MJO in modulating the 20–70-day ISVs in individual summer monsoon regions.
figure 6

Red Arrows indicate lead correlations between MJO convective center in equatorial central Indian Ocean (or negative convective anomaly in tropical western Pacific) and the leading modes of 20–70-day ISVs over individual monsoon regions during a boreal summer and b austral summer. Shading over the ocean denotes the regressed OLR anomaly (W m−2) onto the negative RMM2. The shading over land depicts the leading EOF pattern (mm day−1) of each monsoon (Fig. 3b). The sign of the EOF in Fig. 3b is reversed when the maximum correlation coefficient is negative. The values in red indicate the maximum correlation coefficient, and the values in the bracket show the leading days. All correlations are significant at the 90% confidence level, except for the two African monsoons.

Recent trends

We observe a significant (above 99% confidence level) increasing trend in the amplitude of HFISVs over Southeast China and northern India-Bangladesh during the past four decades (Fig. 7). The largest increasing trend (about 35% per four decades) is over northern India-Bangladesh. In contrast, decreasing trends occur over South America, including significantly (above 95% confidence level) weakening over Venezuela (9%), and Brazil (9%). The LFISVs are enhanced over Central India, Southeast China, and Australia with a trend of 16% (above 90% confidence level), 22% (above 99% confidence level), and 15% (above 90% confidence level), respectively. The decreasing LFISV is found over Brazil (13%, above 99% confidence level).

Fig. 7: Trends in the ISV intensities of individual summer monsoons during 1979–2018.
figure 7

The red bars show the trends (units in percentage change) of the ISV intensity averaged over the wet regions of the leading ISV modes (precipitation anomalies above 0.05 mm day−1 in Fig. 3a) for a 8–20-day ISV and b 20–70-day ISV. The blue bars show the corresponding trends in the area-averaged seasonal-mean precipitation. The ISV intensity is defined by the local-summer standard deviation of filtered precipitation anomaly. Significant trends at the 90% confidence level based on the Mann–Kendall rank statistics are marked by red circles.

The changes in ISV intensity follow those of seasonal-mean precipitation (Fig. 7). The mean precipitation exhibits a significant increasing trend over Southeast China, Central India, northern India-Bangladesh, and Australia. Meanwhile, a significant decreasing trend appears over Venezuela and Brazil, resulting in a significant reduction in ISV intensity. These results suggest that the intensity change of individual monsoon ISVs is determined by local mean-state change. The increased mean precipitation in Asian monsoons during the last four decades is mainly due to the negative phase of the Inter-decadal Pacific Oscillation and the warm phase of the Atlantic Multi-decadal Oscillation, augmented by increased greenhouse gas (GHG) emission63,95.

Implication for sub-seasonal prediction

The extensive database of near real-time ensemble predictions and hindcasts (Table 3), developed for improving prediction skills and physical understanding on the S2S time scale96, provides us an opportunity to investigate sub-seasonal predictability among different land monsoon regions. These prediction ensembles, including 11 S2S prediction models, have been well used to investigate different scientific and modeling issues, as reviewed in de Andrade et al.45. The multi-model mean prediction skill of the 11 S2S prediction systems is focused on in this work (see Methods).

Table 3 Main features of the 11 contributions to the S2S database.

To evaluate the sub-seasonal prediction of the S2S prediction models, the Climate Prediction Center (CPC) precipitation data were used as observations. Although there is a temporal discontinuity at a few grids in the CPC data, the consistency between dominant ISV modes represented by CPC precipitation anomalies (Fig. 3) and those by OLR anomalies (Figs. 4, 5), and between high monsoon ISV-MJO correlations (Fig. 6) and high LFISV contributions over AU, SA, and SAM regions (Table 2), suggests the reliability of using the CPC data to investigate the dominant ISV modes of global monsoons. High-quality observations of monsoon precipitation are necessary to further investigate the characteristics of monsoon precipitation ISV and evaluate model predictability in the future.

For the leading ISV modes over most monsoon regions, the multi-model mean prediction skill of these 11 S2S prediction models decays quickly within the first 3 weeks, from an average correlation of about 0.5 in the first week to below 0.2 in the third (Fig. 8a, b). Divergent prediction skills, however, are observed in individual land monsoon regions. No significant skill can be observed since day 10 for NAF and SAF. EA has very high prediction skills within 1 week, while the skill drops quickly in the second week. AU, however, sustains its significant skill for up to 24 days. The higher prediction skill over Asian and Australian monsoons than the others was also found in the National Centers for Environmental Prediction (NCEP) Climate Forecast System version 297. Among these 11 S2S prediction models, the ECMWF model has the highest prediction skills over most land monsoons (Fig. 8c).

Fig. 8: Divergent sub-seasonal prediction skills over individual summer monsoons.
figure 8

Correlation between observed precipitation anomalies (in CPC data) and multi-model ensemble mean of 11 S2S prediction models’ hindcasts averaged over the wet regions of the leading ISV modes (precipitation anomalies above 0.05 mm day−1 in Fig. 3a) for a NH and b SH summer monsoons. Correlations are calculated from leading 2 to 30 days for the hindcasts initialized around the 1st, 10th, and 20th of each month during the summers of 1999–2010. Shading indicates one standard deviation among the 11 S2S models and dotted lines denote correlation significant at the 90% confidence level. c shows the correlation of weekly precipitation anomalies averaged from Week 2 to Week 4 for multi-model ensemble mean (bar) and each model (symbols).

The multi-model mean prediction skill exhibits high sub-seasonal (Week 2 to Week 4) prediction skills for the leading ISV modes over SA, AU, and SAM (Fig. 8c). The African monsoons, however, have the lowest sub-seasonal prediction skills. This result confirms our hypothesis that the rainfall sub-seasonal predictability of SA, AU, and SAM monsoons with larger LFISV contributions is higher than other monsoons dominated by HFISVs.

Discussion

Our review and investigation show that monsoon regions host the maximum intraseasonal rainfall anomalies during local summer monsoon seasons over the global land areas (Fig. 1a). The averaged ISV intensity (standard deviation) is about half of the seasonal-mean precipitation (9.3 mm day−1) and one-quarter to one-third of the daily precipitation variance (Table 1), and the ISV intensity closely follows the corresponding seasonal-mean precipitation intensity. The area-averaged land monsoon rainfall spectra show a significant 8–20-day (HF) ISV peak in all eight regional monsoons (including the maritime continent) and significant 20–70-day (LF) peaks in AU, SA, and SAM monsoons (Fig. 2). The HFISV accounts for about 53–70% of the total ISV in the land monsoon regions, dominating the sub-seasonal variability of monsoon rainfall, while the LFISV has a relatively large fractional contribution to the total ISV over AU (47%), SA (45%), and SAM (40%) monsoons (Table 2).

The wet phase of ISV in each monsoon region is preceded by diverse precursors and propagation routes, as summarized in Fig. 9, which can be seen as the traces of anomalous anticyclone, southerly, or easterly (marked by letters) in the associated regression maps (Figs. 4, 5 and Supplementary Figs. 4, 5, 7, 8). The Asian monsoon’s HFISVs are primarily related to the westward-propagating quasi-biweekly oscillation disturbances that originate from the western equatorial Pacific, and the EA monsoon is additionally affected by southeastward-propagating mid-latitude wave trains. The NAM and NAF monsoons’ HFISVs find their preceding signals from convectively coupled mixed Rossby-gravity waves and equatorial Kelvin waves, respectively. In contrast, all three SH monsoons are initiated by mid-latitude wave trains. The MJO directly affects the 20–70-day LFISVs of the AU monsoon during its eastward journey in boreal winter and Asian summer monsoon by northward and northeastward propagations. On the other hand, the MJO influences American and African monsoons primarily through its associated Kelvin wave component. Besides, the MJO has a substantial impact on SAM through its mid-latitude teleconnection in the SH. These divergent origins and propagations of different monsoon ISVs are primarily determined by their geographic locations and associated background circulations.

Fig. 9: Schematic diagram of propagation routes for dominant ISV modes of individual summer monsoons.
figure 9

Thin arrows show the propagation routes before the wet peak phase of leading ISV mode (same shading as in Fig. 3) in each monsoon region. The purple and red arrows illustrate the preceding quasi-biweekly oscillation in a and MJO in b, respectively. The eastward and westward blue arrows indicate the preceding equatorial Kelvin and mixed Rossby-gravity waves, respectively. The yellow arrows indicate the mid-latitude intraseasonal Rossby waves. The broad light-blue arrows in a denote the background westerly zones where mid-latitude wave trains prevail.

The HFISV activity has significantly intensified over the past 40 years in Southeast China and northern India-Bangladesh. Similarly, the LFISV displays a significant increasing trend over Central India, Southeast China, and northern Australia (Fig. 7). The top increasing trend reached 35% over northern India-Bangladesh. On the other hand, a significant decreasing trend is observed over Venezuela for the HFISV and over Brazil for both HF and LFISVs. These ISV intensity changes follow the corresponding seasonal-mean precipitation changes, suggesting a local monsoon mean-state control over ISV changes. Observed changes in the mean monsoon rainfall vary by region with significant decadal variations. The NH land monsoon rainfall as a whole has increased since 1980 due to the competing influences of internal climate variability and radiative forcing from GHGs and aerosol forcing; however, it remains a challenge to quantify their relative contributions98.

By comparing HF and LFISVs among different regional land monsoons, our finding of the dominant contribution of HFISV (53–70%) to the total ISV in individual monsoons explains why the prediction skill of the S2S prediction models decays quickly within 3 weeks (Fig. 8a, b), suggesting that better understanding and modeling of HFISV should be critical to sub-seasonal prediction. The new comparison of different origins and propagations associated with HFISVs in each monsoon region may also help identify models’ strengths and weaknesses in reproducing regional characteristics of the ISV: the observed close connections between the NH monsoon HFISVs and the convectively coupled equatorial waves suggest the importance of faithfully simulating convectively coupled equatorial waves in the forecast models; while given the close linkage between the HFISV and mid-latitude wave trains in the SH and EA monsoons, untangling and quantifying the roles of mid-latitude and off-equatorial processes in monsoon rainfall ISV should be a priority.

Our new finding of higher LFISV contributions in SA, AU, and SAM monsoons explains the reason why the S2S prediction models have higher sub-seasonal prediction skills over these three land monsoons than over the others that are dominated by HFISVs (Fig. 8c), which calls for investigation of two-way interaction between regional monsoon and MJO and between regional monsoon ISV and off-equatorial processes. The divergent trends of ISV intensity observed over different land monsoons also challenge our effort to improve sub-seasonal prediction by only focusing on the ISV itself, because the correct simulation of the mean state in a changing climate is also important.

Methods

Observation data

The data used in this study include (1) daily CPC global precipitation data with a high resolution of 0.5° × 0.5°, provided by the National Oceanic and Atmospheric Administration (NOAA)99,100; (2) daily outgoing longwave radiation (OLR) data with a resolution of 2.5° × 2.5° from the NOAA101, as a proxy for large-scale organized deep convection in the tropics12,13,42; (3) daily winds of the Reanalysis II data from the NCEP102; and (4) daily RMM index of Wheeler and Hendon42 from the Australian Bureau of Meteorology. The change of MJO signal over the central equatorial Indian Ocean can be represented by the negative RMM2 index42. Since CPC data were lacking in some grids for a few years, those years are not used in our calculation for these grids.

The study period covers 40 years, from 1979 to 2018. The summer season denotes May to September for the NH and November to March for the SH. Here, we focus on land monsoon precipitation. The global monsoon domains are defined by the regions where the summer-minus-winter precipitation exceeds 300 mm, and summer precipitation exceeds 55% of the annual total103,104.

S2S prediction project database

This study also uses hindcasts of precipitation from 11 global prediction systems participating in the S2S prediction project3. This project was jointly implemented by the World Weather Research Program (WWRP) and World Climate Research Program (WCRP), aiming to improve prediction skills and physical understanding on the S2S time scale105,106. Details of the main features of these S2S model hindcasts are listed in Table 3. To investigate the prediction skill, all CPC and S2S hindcasts are linearly interpolated to a horizontal resolution of 1.5° × 1.5°. Since these S2S models has different hindcast ensemble sizes, an ensemble average is performed for each model thus they have the same weight. Due to different hindcast frequencies, we select three hindcasts from initial conditions around the 1st, 10th, and 20th day of each month for each model. The period of 1999–2010 when all models have hindcasts is studied.

Dominant modes of ISVs

In this work, the seasonal cycle is removed before we calculate the daily anomaly by subtracting the daily climatology. We also apply the Lanczos bandpass filter107 to daily anomalies to extract the HF and LF signals. When the power spectra of precipitation time series are displayed in an area-conserving format, i.e., the logarithm of frequency versus the product of power and frequency, the variance for each frequency band is proportional to the area on this type of spectral map42,108. Thus, the contributions of HF and LFISV signals to the total ISV variance are defined by the ratio of the averaged power spectra of their associated frequency bands versus the total 8–70-day band.

We perform the empirical orthogonal function (EOF) analysis on the filtered data to reveal the dominant patterns of the HF and LFISVs during the summer monsoon season. Each principal component (PC) is normalized by its standard deviation, and the EOF pattern is scaled by multiplying this standard deviation. The Theil-Sen trend estimation method is applied to the linear trend, significantly more accurate than simple linear regression for skew and heteroscedastic data109,110. The statistical significance levels are tested using the Mann–Kendall rank statistics111. Lead-lag regression maps are used to detect potential predictability sources of the ISV patterns for each monsoon region. We test statistical significance based on the two-tailed Student’s t-test, and the significance of the results is assessed at a 90% confidence level without especially mentioned in this study. Since the filtering method is used in this work, the effective degree of freedom is considered to test the significance112, which is defined as \(N_e = N(1 - r1)/(1 + r1)\), where N is the sample size of total days and r1 is the lag-one autocorrelation.