Intraseasonal variability of global land monsoon precipitation and recent change

Accurate prediction of global land monsoon rainfall on a subseasonal (2-8 weeks) time scale has become a worldwide demand. Current forecasts of weekly-mean rainfall in most monsoon regions, however, have limited skills beyond two weeks. Given that two-thirds of the world’s population lives in the monsoon regions, this challenge calls for a more profound understanding of monsoon intraseasonal variability (ISVs). Our comparison of individual land monsoons shows that the high-frequency (HF; 8-20 days) ISV, crucial for the Week 2 and Week 3 predictions, accounts for about 53-70% of the total (8-70 days) ISV in various monsoons, and the low-frequency (LF; 20-70 days) ISV has a relatively high contribution over Australia (AU; 47%), South Asia (SA; 43%), and South America (SAM; 40%) monsoons. The leading modes of HFISVs in Northern Hemisphere (NH) monsoons primarily originate from convectively coupled equatorial Rossby waves (Asia), mixed Rossby-gravity waves (North America, NAM), and Kelvin waves (northern Africa, NAF), while from mid-latitude wave trains for Southern Hemisphere (SH) monsoons and East Asian (EA) monsoon. The Madden-Julian Oscillation (MJO) directly regulates LFISVs in the Asian-Australian monsoon while affecting the American and African monsoons by exciting Kelvin waves and mid-latitude teleconnections. During the past four decades, the HF (LF) ISVs have considerably intensied over the Asian (Asian-Australian) monsoon but weakened over the American (SAM) monsoon. Subseasonal-to-seasonal (S2S) prediction models do exhibit higher subseasonal (Weekly 2-Weekly 4) prediction skills over SA, AU, and SAM monsoons that have larger LFISV contributions than the other monsoons. The results suggest an urgent need to improve the simulation of convectively coupled equatorial waves and two-way interactions between regional monsoon ISVs and mid-latitude processes and between MJO and regional monsoons, especially under the global warming scenarios.


Introduction
High-quality social development with effective disaster prevention, logistics-planning, agriculture production, and decision-making requires accurate subseasonal rainfall prediction over land monsoon regions where about two-thirds of the world's population live 1 . Tropical intraseasonal variabilities (ISVs) and the associated teleconnection are dominant sources of subseasonal predictability 2-4 .
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 subseasonal prediction skill and predictability estimates has primarily focused on the leading modes of the MJO or boreal summer intraseasonal oscillation 3,[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) index 42,43 . However, for up to four weeks, the MJO's rainfall prediction skill is con ned to a narrow equatorial belt, and the rainfall prediction skills drop rapidly beyond two weeks in vast off-equatorial monsoon regions 44,45 , as shown in Fig. 1. User community concerns with weekly rainfall prediction rather than the MJO index. What limits the subseasonal predictability of monsoon rainfall? This is the forefront challenge faced by the subseasonal prediction community.
Why does the weekly-mean monsoon rainfall prediction skill decline 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 variability 46 . 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 the HFISV in subseasonal prediction and predictability. The HFISVs in various monsoon regions have been documented. Over SA, HFISV was dominated by the westward propagating quasi-biweekly oscillation 16,[47][48][49][50] . The signi cant HFISV of EA precipitation was linked to preceding midlatitude wave trains and tropical quasi-biweekly oscillation 18,[50][51][52] . HFISVs were also reported over NAM 53 , NAF 23,24,54 , AU 27,55 , and SAM 29,33,[56][57][58] . 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 con guration and involves different atmosphereocean-land interaction processes and convective and mesoscale systems 59 . 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 eastwardpropagating planetary-scale MJO circulation system and mid-latitude disturbances may signi cantly regulate and coordinate the regional ISVs by changing monsoon circulations. These coordination and regional in uences 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 subseasonal prediction of regional land monsoon rainfall.
Our particular interests also include examining recent changes in regional monsoon ISVs under global warming. The seasonal-mean rainfall in the NH monsoon regions has increased considerably since 1979 60 , and the MJO residence time was observed to increase over the Indo-Paci c Maritime Continent by 5-6 days during the last four decades due to the warm pool expansion 61 . These changes of mean state and MJO might have in uenced the land monsoon ISV. The regional ISV was enhanced since 1993 over South China 62 . However, the changes of both HF and LF ISV intensities in most land monsoon regions 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 subseasonal prediction skills over different land monsoon regions in current subseasonal-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 subseasonal prediction of land monsoon rainfall.

Results
Relative importance of HF and LF ISVs in summer monsoon regions The maximum intraseasonal rainfall anomalies tend to anchor over the global land monsoon regions (Fig. 2a). The average ISV intensity of the global land monsoon rainfall, measured by the standard deviation of the 8-70-day ltered 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 onequarter 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.
The regionally averaged land monsoon precipitation spectra show a signi cant 8-15-day HFISV peak in NH monsoons and a signi cant 10-20-day HFISV in SH monsoons (Fig. 3). The precipitation spectra also show signi cant LF peaks in AU (30-60 days) and in SA and SAM (20-40 days). No signi cant signal is observed beyond 70 days. Thus, we de ne the HF component with a time scale of 8-20 days (~2-3 weeks) and the LF component, 20-70 days (~4-9 weeks).
The relative contributions of these HF and LF ISVs to the total ISV variance are de ned 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 Australian 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, where MJO has less in uence, the HFISV accounts for about 66%-70% of the total ISV variance, dominating the ISV. The results suggest that MJO has a more signi cant in uence 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 AU monsoons have 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 amplitude 34,63,64 , 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 nd out whether the rainfall subseasonal predictability in these monsoon regions is also signi cantly lower than in AU, SA, and SAM.
The Maritime Continent is part of the monsoon region in terms of its annual reversal of the prevailing wind 65,66 . The seasonal-mean precipitation and 8-70-day ISV intensity averaged over the Maritime Continent (10 o S-5 o N; 100 o -150 o E) are 12.5 mm day -1 and 5.8 mm day -1 , respectively, during its wet season from November to March ( Supplementary Fig. 3a). During the wet season, the Maritime Continent rainfall only shows a signi cant 8-20-day HFISV component (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 rst identify the leading spatial pattern and activity center of the HFISV of summer precipitation in each monsoon region. Figure 4a 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 in the tropical monsoon trough and the other in the subtropical convergence zone (poleward of 20 o latitude). The uniform precipitation anomalies in the NAF, AU, SAM, and the 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 rst ones (Supplementary Table 1). Together, the two leading modes explain about 15% of the variance.
Where are the peak phases of the leading 8-20-day ISV modes of land monsoon summer precipitation originated? Following many previous works 49,52,67 , we use the sequential leading maps of regressed OLR anomalies associated with these leading ISV modes 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 (Fig. 5). On Day 0, the regressed OLR anomalies are consistent with the precipitation anomalies associated with the leading 8-20-day ISV pattern shown in Fig. 4a.
Over SA (Fig. 5a), the low-level anticyclonic anomaly that couples with the suppressed convection This feature is consistent with previous ndings 16,47,48 . An apparent upper-tropospheric wave train also moves westward along 40 o N 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 signi cantly modulates the mid-latitude circumglobal teleconnection 68 , "silk road" teleconnection 69,70 , and the Tibetan Plateau 71 on the HF intraseasonal time scale.
In EA (Fig. 5b) 10 . Note that from Day -2 to Day 0 the convection over Southeast China is enhanced by the convergence of southwesterly wind anomalies to its northwest, forming the dipolar structure. In addition to the tropical origin mentioned above, a lowertropospheric cyclonic wind anomaly emerges to the east of Lake Baikal and propagates southeastward ( Fig. 5b). 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 Paci c. The evidence here seems to con rm previous ndings concerning the in uence of the mid-latitude wave trains on the quasi-biweekly variability of EA summer monsoon 51,52,74 .
Over the NAM and NAF monsoon regions, the OLR anomalies are relatively weak (Figs. 5c and 5d). However, their precipitation anomalies are not (Fig. 4a), suggesting that shallower convection may prevail over these two regions, a feature similar to the corresponding seasonal-mean precipitation [75][76][77] . 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. 5c). 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 ow looks like a mixed Rossby-gravity wave 78 . We suggest that the 8-20-day ISV in the NAM region originates from an equatorial mixed Rossbygravity wave and gradually transforms to a convectively coupled Rossby wave. No signi cant midlatitude 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. 5d). The upper-level Kelvin wave easterly anomalies emerge over the eastern equatorial Paci c around 120 o W on Day -6, propagate eastward, and reach the eastern Atlantic around 10 o W on Day 0, with a phase speed of about 24 m s -1 . The convection over NAF is then enhanced by the divergent ows 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 signi cant mid-latitude wave trains in the lower troposphere (Fig. 6) 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 was originated in the southeastern Indian Ocean on Day -8 and propagated eastward to Southeast Australia on Day 0 (Fig. 6a). In SAM, a cyclonic anomaly emerges over the South Atlantic along 50 o S on Day -8, which propagates northward rst, and then turns westward on Day -6, triggering the convection over SAM on Day -2 (Fig. 6b). The convection anomaly couples with the cyclonic anomaly and matures on Day 0. The northward and westward propagations were also associated with the sub-monthly oscillation of the Brazil rainfall or South Atlantic Convergence Zone 56,57,79 . The role of mid-latitude wave trains on the ISV is also observed for the SAF monsoon (Fig. 6c). 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. 6 and Supplementary Fig. 5). However, once coupled with convection in the monsoon region, the barotropic structure transforms to a baroclinic structure.
In summary, the ISVs in Asian monsoons are related to the preceding tropical quasi-biweekly oscillation originated from the western equatorial Paci c. The EA monsoon is additionally associated with southward-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 Rossbygravity 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. 52 . 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-20day) ISVs for each monsoon (Fig. 4). Nevertheless, minor differences exist. In SA, the Central Indian precipitation anomaly is stronger than that in northern India-Bangladesh (Fig. 4b). In NAM, the Mexican precipitation anomaly is stronger than the Venezuela anomaly. The similarity between the spatial patterns of the HF and LF ISVs implies that the mean regional monsoon circulation might have similar controls on the HF and LF ISVs. The HF and LF ISVs may also interact with each other.
In the SA, AU, and SAM regions, where the LFISVs have signi cant contributions to the total ISV (Table 2), the rst 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 be geographically-locked ampli cation. This geographically-locked ampli cation of the ISV has also been observed before for SA 80 and for SAM 81,82 , which means that the mean regional monsoon circulation has signi cant control on the LFISVs. Conversely, in other monsoon regions where the HF dominates the ISV, the rst two modes of the LFISVs have comparable fractional variances, suggesting a propagating mode. It would be interesting to explore in the future, at what time scale, the mean monsoon circulation and tropical or mid-latitude intraseasonal forcing have signi cant forcing on the monsoon ISV.
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 Paci c suppressed ISO that originates from the Indian Ocean ( Supplementary Fig. 7b), consistent with previous works 8,48,83 . 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 28,84 . The MJO modulates other monsoon regions primarily by exciting convectively coupled Kelvin waves. Signi cant 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).
Signi cant preceding mid-latitude wave trains are observed in SAM. The wave train associated with SAM originated to the east of Australia, tied to the MJO convection anomaly over the tropical western Paci c ( Supplementary Fig. 8c). This Indo-Paci c MJO's impact on SAM through the mid-latitude teleconnection has also been observed previously 34,37,56,57,64,[85][86][87] . Since the group speed of the teleconnection is quite fast, and the MJO's signal can reach the America in about one week 88,89 , so that the SAM only lags the Indo-Paci c 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 works 90, 91 .
In summary, the 20-70-day ISVs of the Asian-Australian monsoon are directly affected by the MJO propagation. The Indo-Paci c MJO in uences the American and African monsoons through its eastwardpropagating Kelvin component and also has robust impacts on the SAM through its mid-latitude teleconnection in the SH. Figure 7 summarizes the relationship between the MJO peak wet phase over the central equatorial Indian Ocean (or dry phase over the tropical western Paci c) and each region's leading mode of 20-70-day ISV. All these correlations are signi cant, except for the 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 Paci c, 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. 2b) and the higher contribution of the LFISV to the total ISV (Table 2). In addition, the enhanced MJO convection over the central equatorial Indian Ocean tends to lead the wet Southeast Asia and dry South China by 22 days. The dry anomaly over the tropical western Paci c 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 subseasonal prediction of each regional monsoon.

Recent trends
We observe a signi cant (p<0.01) increasing trend in the amplitude of HFISVs over Southeast China and northern India-Bangladesh during the past four decades (Fig. 8). The largest increasing trend (about 35% per four decades) is over northern India-Bangladesh. In contrast, decreasing trends occur over South America, including signi cantly (p<0.05) weakening over Venezuela (9%) and Brazil (9%). The LFISVs are enhanced over Central Indian, Southeast China, and Australia with a trend of 16% (p<0.1), 22% (p<0.01), and 15% (p<0.1), respectively. The decreasing LFISV is found over Brazil (13%, p<0.01).
The changes of ISV intensity follow those of seasonal-mean precipitation (Fig. 8). The mean precipitation exhibits a signi cant increasing trend over Southeast China, Central India, northern India-Bangladesh, and Australia. Meanwhile, a signi cant decreasing trend appears over Venezuela and Brazil, resulting in a signi cant reduction in ISV intensity. The increased mean precipitation in Asian monsoons during the last four decades is mainly due to the negative phase of the Inter-decadal Paci c Oscillation and the warm phase of the Atlantic Multi-decadal Oscillation, augmented by increased greenhouse gas (GHG) emission 60,92 . The year-to-year variations of the ISV intensity over different land monsoons are also signi cantly correlated to the variations of corresponding mean precipitation (Supplementary Fig. 9). All these results indicate that the intensity change of individual monsoon ISVs is determined by local meanstate change.

Implication for subseasonal prediction
The extensive database of near real-time ensemble predictions and hindcasts (Table 3), developed for improving prediction skill and physical understanding on the S2S time scale 93 , provides us an opportunity to investigate subseasonal predictability among different land monsoon regions. These prediction ensembles, including 11 S2S prediction models, have been well used to investigate different scienti c and modelling issues, as reviewed in de Andrade et al. 45 . The multi-model mean of the 11 S2S prediction systems is focused in this work (see Methods).
To evaluate the subseasonal predictability of the S2S prediction models, the Climate Prediction Center (CPC) precipitation data are used as observations. Although there is temporal discontinuity at a few grids in the CPC data, the consistence between dominant ISV modes represented by CPC precipitation anomalies (Fig. 4) and those by OLR anomalies (Figs. 5 and 6), and between high monsoon ISV-MJO correlations (Fig. 7) 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 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 rst three weeks, from an averaged correlation of about 0.5 in the rst week to below 0.2 in the third (Figs. 9a and 9b). Divergent prediction skills, however, are observed in individual land monsoon regions. No signi cant skill can be observed since day 8 for NAF and since day 10 for SAF (Fig. 9a). EA has very high prediction skill within one week, while the skill drops quickly in the second week. AU, however, sustains its signi cant skill 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 2 94 . Among these 11 S2S prediction models, the ECMWF model has the highest prediction skills over most land monsoons (Fig. 9c). The multi-model mean of the 11 S2S prediction models, thus, exhibits high subseasonal (Week 2 to Week 4) prediction skills for the leading ISV modes over SA, AU, and SAM (Fig. 9c). The African monsoons, however, have the lowest subseasonal prediction skills. This result con rms our hypothesis that the rainfall subseasonal predictability in SA, AU, and SAM monsoons with larger LFISV contributions is higher than the other monsoons dominated by HFISVs.

Summary
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. 2a). The averaged ISV intensity (standard deviation) is about half of the seasonal-mean precipitation (9.3 mm day -1 ) and onequarter 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 signi cant 8-20-day (HF) ISV peak in all eight regional monsoons (including the maritime continent) and signi cant 20-70-day (LF) peaks in AU, SA, and SAM monsoons (Fig. 3). The HFISV accounts for about 53-70% of the total ISV in the land monsoon regions, dominating the subseasonal 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. 10. The Asian monsoon's HFISVs are primarily related to the westwardpropagating quasi-biweekly oscillation disturbances originated from the western equatorial Paci c, and the EA monsoon is additionally affected by southward-propagating mid-latitude wave trains. The NAM and NAF monsoons' HFISVs nd their preceding signals from convectively coupled mixed Rossby-gravity waves and equatorial Kelvin waves, respectively. In contrast, all three SH monsoons are initiated by midlatitude 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 in uences American and African monsoons primarily through its associated Kelvin wave component. Besides, the MJO has a substantial impact on the SAM through its mid-latitude teleconnection in the SH. These divergent origins and propagations associated with different monsoon ISVs are primarily determined by the geographic locations of the monsoons and their associated background circulations.
The HFISV activity has signi cantly intensi ed over the past 40 years in Southeast China and northern India-Bangladesh. Similarly, the LFISV displays a signi cant increasing trend over Central India, Southeast China, and northern Australia (Fig. 8). The top increasing trend reached 35% over northern India-Bangladesh. On the other hand, a signi cant decreasing trend is observed over Venezuela for the HFISV and over Brazil for both HF and LF ISVs. 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 signi cant decadal variations. The NH land monsoon rainfall as a whole has increased since 1980 due to the competing in uences of internal climate variability and radiative forcing from GHGs and aerosol forcing; however, it remains a challenge to quantify their relative contributions 95 .
By comparing HF and LF ISVs among different regional land monsoons, our nding of 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 three weeks (Figs. 9a and 9b), suggesting that better understanding and modeling of HFISV should be critical to subseasonal 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 midlatitude 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 nding of higher LFISV contributions in SA, AU, and SAM monsoons explains why the S2S prediction models have higher subseasonal prediction skills over these three land monsoons than over the others that are dominated by HFISVs (Fig. 9c), 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 subseasonal prediction by only focusing on the ISV itself, because correct simulation of the mean state in a changing climate is also important.

Methods
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) 96,97 ; 2) daily outgoing longwave radiation (OLR) data with a resolution of 2.5° × 2.5° from the NOAA 98 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, 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 1 st , 10 th , and 20 th day of each month for each model. The period of 1999-2010 when all models have hindcasts is studied.
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 lter 104 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 map 42,105 . Thus, the contributions of HF and LF ISV signals to the total ISV variance are de ned by the ratio of the averaged power spectra of their associated frequency bands versus the total 8-70-day ISV band.
We perform the empirical orthogonal function (EOF) analysis on the ltered data to reveal the dominant patterns of the HF and LF ISVs 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, signi cantly more accurate than simple linear regression for skew and heteroscedastic data 106,107 . The statistical signi cance levels are tested using the Mann-Kendall rank statistics 108 . Lead-lag regression maps are used to detect potential predictability sources of the ISV patterns for each monsoon region. We test statistical signi cance based on the two-tailed Student's t-test. Since ltering method is used in this work, the effective degree of freedom is considered to test the signi cance 109 , which is de ned as where N is the sample size of total days and r1 is the lag-one autocorrelation.

DATA AVAILABILITY
The data that support the ndings of this study are openly available and they can be found in the respective references. Tables   Table 1. The averaged ISV intensity for each regional land summer monsoon precipitation measured by the standard deviation shown in Fig. 1. Also shown is the fractional variance contribution of ISV to the total daily variance.    Power spectra of ISVs of individual monsoons. Power spectra of the 3-day running mean precipitation anomalies (mm day-1) averaged over each monsoon region for 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% con dence levels.  shown. Letter A tracks the center of the anomalous anticyclone; C, the cyclone; S, the southerly; N, the northerly; E, the easterly; and W, the westerly.  onto the negative RMM2 and the leading EOF pattern (mm day-1) of each monsoon shown in Fig. 4b, respectively. The sign of the EOF in Fig. 4b is reversed when the maximum correlation coe cient is negative. The values in red indicate the maximum correlation coe cient, and the values in the bracket show the leading days, respectively. All correlations are signi cant at the 90% con dence level, except for the two African monsoons. The blue bars show the corresponding trends in the area-averaged seasonal-mean precipitation. The ISV amplitude is de ned by the local-summer standard deviation of ltered precipitation anomaly. Signi cant trends at the 90% con dence level based on the Mann-Kendall rank statistics are marked by red circles.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. Supplementary.docx