Research ArticleCLIMATOLOGY

Skillful prediction of summer rainfall in the Tibetan Plateau on multiyear time scales

See allHide authors and affiliations

Science Advances  09 Jun 2021:
Vol. 7, no. 24, eabf9395
DOI: 10.1126/sciadv.abf9395


Skillful near-term climate predictions of rainfall over the Tibetan Plateau (TP), the Asian water tower, benefit billions of people. On the basis of the state-of-the-art decadal prediction models, we showed evidence that although the raw model outputs show low predicted ability for the summer Inner TP (ITP) rainfall due to low signal-to-noise ratios in models, we can produce realistic predictions by extracting the predictable signal from large ensemble predictions along with a postprocessing procedure of variance adjustment. The results indicate that the summer ITP rainfall is highly predictable on multiyear time scales. The predictability of ITP rainfall originates from the Silk Road pattern driven by sea surface temperature over the subpolar gyre region in North Atlantic. Real-time forecasts suggest that the ITP will become wetter, with 12.8% increase in rainfall during 2020–2027 relative to 1986–2005. Our results will help the water resources management in the surrounding regions.


The Tibetan Plateau (TP), known as the “roof of the world” or the “3rd pole of the earth” (1, 2), has the most glaciers outside the Arctic and Antarctic (3). The meltwater feeds more than 10 prominent rivers, including the Yangtze River, the Yellow River, and the Ganges River, supplying water resource to more than 2.0 billion people, playing a role of “the Asian Water Tower” (Fig. 1A) (4). The enormous latent heating released from the precipitation along with the surface sensible heating makes the TP a heat source in boreal summer [June-July-August (JJA)], subsequently influencing regional and even global climate, via interactions with the atmospheric circulations (5). A rapid warming has been witnessed in the TP in recent decades along with retreated glaciers and degraded permafrost (6), which thereby threatens water supplies to the surrounding regions (7). Hence, forecasting future climate change is imperative for the adaptation and mitigation planning in this ecologically vulnerable region.

Fig. 1 The TP and the CMIP6 multimodel ensemble projections of the TP summer precipitation.

(A) The surface elevation of the TP (shaded; unit: m). Distribution of natural lakes (blue shades) and major rivers (light blue lines; including Yellow, Yangtze, Mekong, Salween, Brahmaputra, Ganges, Indus, and Amu Darya) on the TP are also shown. The TP is divided into 12 large river basins (26). (B) Time series of observations and CMIP6 multimodel ensemble projections of the summer precipitation anomalies (unit: mm/day) over the TP relative to 1986–2005. The historical simulations (gray lines) during 1960–2014 and projections under SSP1-2.6 (green lines), SSP2-4.5 (blue lines), SSP3-7.0 (purple lines), and SSP5-8.5 (red lines) emission scenarios during 2015–2100 from 19 CMIP6 models (table S1) are used. The thick (dashed) lines represent the ensemble mean (individual ensemble) simulations, respectively. The observations were derived from Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE) (black line), Global Precipitation Climatology Project (GPCP) (orange line), and Global Precipitation Climatology Centre (GPCC) (brown line). Both the observations and simulations are smoothed with a 9-year running average. The inset figure in (B) emphasizes the projections for the period of 2014–2040.

The urgent need for the timely and effective near-term (2021–2040) climate information has been recognized by the government and scientists in recent decades (8, 9). Previous studies have provided the multidecadal to century-scale climate change projections based on the coupled general circulation models (CGCMs) from the Coupled Model Intercomparison Project (CMIP) along different pathways, which take account of the response to external radiative forcing, including the anthropogenic and natural factors (10). The CGCMs from CMIP5 have projected an increased temperature over TP by about 1.4° and 1.6°C in 2006–2050 under the Representative Concentration Pathway (RCP) RCP4.5 and RCP8.5 emission scenarios, respectively (11). The greatest warming is projected in winter, and the warming rates are amplified with elevation (12). Compared with the temperature, there is much uncertainty in the projections of the TP precipitation during the rainy season (JJA), changing with a range from −1.8 to 15.2% in 2016–2035 under the RCP4.5 (13). Nevertheless, beyond the long-term climate changes determined by external radiative forcing, near-term climate changes are affected by natural internal climate variability (14). For example, compared to the long-term period (2081–2100), the near-term projection of TP summer rainfall under different emission scenarios is indistinguishable from each other because of internal variability (Fig. 1B). So far, how to narrow down the uncertainty in the near-term climate change projection remains an open question (15).

The near-term climate prediction (16) is a valuable tool to overcome the defect of projection uncertainties. The near-term climate prediction, also referred to as decadal climate prediction, is an emerging research field in climate science in recent decades (17). Skillful decadal predictions can bridge the gap between the seasonal-to-interannual predictions and long-term projections (18). Decision-makers concerned with the climate adaptation and resilience could benefit greatly from the skillful decadal climate prediction. In general, both externally forced and internally generated climate variability contribute to decadal climate predictability (19). Decadal climate predictions can be enabled by initialized climate model simulations, which reproduce the historical externally forced variability by prescribing the external radiative forcing as in the historical climate simulations and forecast the internally generated components of the climate system by using data assimilation systems to preset the observed climate state at the beginning of the prediction. Initialized predictions have shown that the global mean surface air temperature (20) and the sea surface temperatures (SSTs) in the North Atlantic Ocean (21), the Indian Ocean (22), and the Southern Ocean (23) can be well predicted in the hindcasts or the retrospective forecast, while the confidence in precipitation predictions over most land regions is generally lower (24).

Precipitation changes in the TP can cause severe geological hazards (25). For instance, the Inner TP (ITP; also referred to as Qiang-tang Plateau), with the largest number of natural lakes, has experienced drastic lake expansion since the mid-1990s (26), in conjunction with a notable decadal variations of boreal summer precipitations (27). How the boreal summer ITP rainfall will change in the near future remains unknown. Here, on the basis of the decadal hindcast experiments from 10 state-of-the-art decadal prediction models of the CMIP6 Decadal Climate Prediction Project (DCPP) (28), we showed evidence that the summer ITP rainfall is highly predictable on multiyear time scales. The predictable signal is extracted from large ensemble predictions along with a postprocessing procedure of variance adjustment. Real-time forecasts indicate that the ITP will experience a 12.8% increase in summer mean precipitation during 2020–2027 relative to 1986–2005.


Decadal prediction skill

To evaluate the decadal prediction skill for boreal summer TP rainfall, we use the Global Precipitation Climatology Centre (GPCC) monthly precipitation dataset as the observation, which is consistent with other frequently-used rainfall datasets with respect to the interdecadal variability (figs. S2 and S3). We calculate the anomaly correlation coefficient (ACC) skill score and the mean squared skill score (MSSS) (see Materials and Methods) for boreal summer rainfall for the predictions averaged over 2 to 9 years in DCPP multimodel ensemble predictions. In total, 138 ensemble members from 10 decadal prediction models are used (table S2). We find significant ACC skills in the main body of the TP by using the 138-member ensemble mean prediction (Fig. 2A). The ACC values are greater than 0.7 over the central TP, with the center located in the ITP. Hence, the ensemble mean prediction is able to capture the phase of decadal variability of the boreal summer ITP rainfall. We further estimate the relative contributions of the initialized and uninitialized components to the overall ACC skill in the ensemble mean prediction (see Materials and Methods and fig. S4). For the average over the ITP, 86.7% (13.3%) of total skill is attributed to the uninitialized (initialized) components of the prediction, respectively (see Materials and Methods), implying that the external radiative forcing has great influence on the predictions.

Fig. 2 Decadal prediction skill for 2- to 9-year forecast periods of the averaged summer rainfall over the TP.

(A) ACC skill for 138-member ensemble mean prediction. (B) ACC skill for 10-member ensemble mean prediction (averaged over all possible combinations). (C) MSSS skill for 138-member ensemble mean prediction. Only the values higher than 0 are shown. (D) RPC calculated from 138-member predictions. The dots denote values passing the 95% confidence level (see Materials and Methods). The elevation of 2500 m is used as the outline of TP. The outline of the ITP is also shown.

Although high ACC skill is seen in 138-member ensemble mean prediction, the amplitude of the 138-member mean predictions is small, indicated by the low MSSS skill (Fig. 2C). The ACC skill calculated by the 10-member ensemble mean is lower (Fig. 2B), suggesting that a small ensemble has a limited ability to predict the phase of decadal variability. Both the evolutions of the actual climate system and the model prediction are composed of predictable (the signal) and unpredictable (the noise) components. The inconsistence between ACC and MSSS in the ensemble mean predictions and the requirement for a large ensemble to obtain higher ACC arise from the low signal-to-noise ratio in models compared with the observations. The disagreements in the signal-to-noise ratio between models and observations can be equivalently measured by the ratio of predictable components (RPCs) (see Materials and Methods). Ideally, the RPC should be one for perfect models. We find a higher-than-expected RPC over the central TP, with values significantly greater than one (Fig. 2D). High values of RPC indicate that the models with low signal-to-noise ratio underestimate the predictability of the real world (29). The ensemble members contain much unpredictable noise and weak predictable signal (30, 31). Therefore, the ensemble mean predictions conducted by larger ensemble can effectively eliminate the noise component and extract the weak predictable signal. The ACC is not affected by the magnitude of the predictable signal, and the significant ACC arises from the ensemble mean predictions. In contrast, the MSSS skill, which measures the amplitude of variability, is low because of the weak predictable signal. It is also evident that the spatial distribution of RPC is consistent with the regions that have added value of ACC skill with increased ensemble sizes (Fig. 2, A and B), further demonstrating that large ensemble is needed to obtain a skillful prediction.

In summary, on the basis of the decadal hindcast experiments from 10 state-of-the-art decadal prediction models of the CMIP6 DCPP, we show evidence that the models have lower signal-to-noise ratio compared with the observations for the prediction of the boreal summer rainfall over the central TP. As a result, the 138-member ensemble mean has a high ACC skill but a low MSSS skill, and the ACC skill increases with the ensemble size.

Prediction of ITP rainfall and variance adjustment

We focus on the prediction of summer rainfall in the ITP, where there exists the largest number of natural lakes over TP (32). We see skillful predictions in the ITP region (Fig. 2A). In the observations, the summer rainfall over ITP shows strong decadal variability, with a decrease from the 1960s to the 1990s and a rapid increase thereafter (Fig. 3A, black line).

Fig. 3 Decadal prediction for the summer ITP precipitation.

(A) The 2- to 9-year ensemble mean predictions of the summer ITP precipitation (unit: mm/day) from decadal hindcasts (raw hindcast ensembles, red line), the quasi–real-time predictions (raw prediction ensembles, blue line), and the corresponding observations (black line). (B) Same as in (A), but after applying variance adjustment to the raw predictions. Cor noted in parentheses is the correlation coefficient between the observation and the predictions. Shading in (A and B) show the 5 to 95% uncertainty range. (C) Correlation of the predicted summer ITP rainfall with the observational reference along the forecast time for 4-year averages. The black dashed line represents the persistence predictions. The vertical segments denote 5 to 95% uncertainty range of all possible combinations. (D) Correlation as a function of ensemble size for model ensemble mean predictions of an independent ensemble member (blue line; averaged over all possible combinations) and for model ensemble mean predictions of the observations (red line). Shading shows the 5 to 95% confidence interval. Time series of the predictions and the observations in (C and D) are detrended. (E) Relationship between the 8-year running mean of the summer ITP rainfall and the summer subpolar gyre (SPG) SSTA in observations. (F) Same as in (E), but for the predicted 2- to 9-year summer SPG SSTA and the summer ITP rainfall from the decadal hindcasts. The COR noted in the top left in (E and F) is the correlation coefficient for the observations and ensemble mean predictions, respectively.

We compare the time series of the ITP rainfall index derived from the raw model outputs to the observations (Fig. 3A, red line). The raw 138-member ensemble mean prediction has a high correlation coefficient with the observation (ACC = 0.74, significant at the 1% level), but the predicted magnitude is weaker than the observation because of the low signal-to-noise ratio in the models (RPC = 4.05). The total variability of the 8-year averaged ITP summer precipitation in individual model members measured by SD is 0.07 to 0.22 mm/day (5 to 95% range and 2- to 9-year prediction), which is not significantly different from the observations (0.14 mm/day). Hence, the variance of the predictable signal in the observations is about 15 times larger than that in the models. In this case, the models can provide better predictions for the real world than for the model itself (30). In the model’s predictions for the real world, the prediction skills of the summer ITP rainfall gradually increase with the increasing ensemble size, with ACC reaching the saturation point at 0.77 (Fig. 3D, red line), while the skills for the model in predicting one realization of the model itself are lower, and the saturation point of ACC is 0.19 (Fig. 3D, blue line). This contrasts the “perfect model” assumption that low skill is expected from the low signal-to-noise ratio in models and thus is referred to as the “signal-to-noise ratio” paradox (31). By using large ensemble, we find a steady increase of skill along with increasing ensemble sizes for the forecast periods ranging from 1 to 4 years to 7 to 10 years (Fig. 3C).

To improve the prediction, we need to rescale the variance of the predictable signal and the unpredictable noise in the raw predictions. Here, we use the observations to adjust the raw predictions with a proper consideration of the influence of long-term linear trend (see Materials and Methods). After variance adjustment, the predictions are highly consistent with the observations, as evidenced by the ACC value increasing to 0.85 (significant at the 1% level). The prediction reasonably captures the decrease from the 1960s to the 1990s and the rapid increase in the late 1990s of the summer ITP rainfall (Fig. 3B).

In summary, we find that the signal-to-noise ratio for the prediction of the summer ITP rainfall in models is low. The raw model outputs have a high ACC skill but a weak predicted magnitude. To produce realistic predictions, we extract the predictable signal using a large ensemble and apply a postprocessing procedure of variance adjustment. The results indicate that the summer ITP rainfall is highly predictable on multiyear time scales.

Sources of decadal predictability

What are the sources of TP rainfall decadal predictions? We examine the relationship between the 8-year running averaged summer SST anomalies (SSTAs) of the subpolar gyre (SPG) region in the North Atlantic (50° to 65°N; 60° to 10°W) (33) and the summer ITP rainfall in both the observation and the prediction (Fig. 3, E and F). In the observations, the SPG SSTA has a close relationship with the summer ITP rainfall, with the correlation coefficient reaching 0.87 (significant at the 1% level). In the predictions, the correlation coefficients between the predicted SPG SSTA and the predicted summer ITP rainfall range from −0.265 to 0.545 among the ensemble members, while the correlation in the ensemble mean prediction is 0.880 (significant at the 1% level). A warmer (colder) SPG SSTA is favorable for increased (decreased) summer rainfall over the ITP on the interdecadal time scale both for the observations and decadal hindcasts. Hence, the decadal predictability of ITP rainfall originated from the North Atlantic.

To reveal the physical processes underpinning the decadal predictability of boreal summer ITP rainfall, we investigate the SSTA and atmospheric circulation anomalies associated with the variations of summer ITP rainfall in the upstream regions. The simultaneous regressions of 200-hPa eddy geopotential height (see Materials and Methods) show a zonal oriented wave structure over the mid-latitude Eurasian continent, and the wave activity fluxes (WAFs) for stationary Rossby waves propagate from the North Atlantic Ocean to the downstream East Asian regions along the upper-level jet stream waveguides (Fig. 4A). We see a significant warming over the SPG region (Fig. 4B), indicating that the decadal change of ITP summer rainfall is linked to the Atlantic multidecadal variabilities (34).

Fig. 4 Mechanism for the decadal variability of ITP summer precipitation and sources of decadal predictability in DCPP.

Regressions onto the 8-year running averaged summer ITP rainfall index for (A) the observed 200-hPa eddy geopotential height (shaded; unit: m) and WAF (vector; unit: m2·s−2) and (B) the observed SSTA (shaded; unit: °C) and land precipitation (shaded; unit: mm/day). (C and D) Same as in (A and B), but for the regressions of the 2- to 9-year averaged ensemble mean predictions onto the 2- to 9-year averaged ensemble mean predicted ITP summer precipitation index. The dots denote values passing the 95% confidence level (see Materials and Methods). (E and F). ACC skill for 2- to 9-year forecast periods of the averaged (E) 200-hPa eddy geopotential height and (F) SSTA and land precipitation. The dots denote values passing the 95% confidence level (see Materials and Methods).

To further demonstrate the key role of the SPG SSTA in driving the wave trains, we use the maximum covariance analysis (MCA) to extract the maximum coordinated variation of SSTA over the SPG region and the summer 200-hPa meridional wind anomalies over the mid-latitude regions over the Eurasian continent (see Materials and Methods). The leading MCA mode is highly consistent with the circulation anomalies related to the summer ITP rainfall, indicating that the SPG SSTA and the wave trains over the Eurasian are tightly coupled on the decadal time scale (fig. S5). The wave pattern over the Eurasian is the interdecadal Silk Road pattern (3537), which is a typical teleconnection pattern in the summertime that encompasses several geographically fixed centers along the upper tropospheric subtropical jet over the Eurasian continent.

How does the interdecadal Silk Road pattern influence the summer ITP rainfall? Two possible mechanisms are suggested. The first is related to the Silk Road pattern–related anomalous upper troposphere anticyclone over northeast of the TP under the influence of the positive phase of the Atlantic multidecadal oscillation (38). The anomalous anticyclone can weaken the westerly winds in the lower levels, prevent the export of water vapor to the eastern boundary of the TP, and thus affect the summer ITP precipitation variability. The second mechanism emphasizes the role of the Silk Road pattern–related upper troposphere anomalous cyclone to the west of the TP (27), which facilitates water vapor meridionally intruding from the Arabian Sea into the ITP and thus also influences the variations of the summer ITP rainfall. Hence, both the upper troposphere anticyclone over the northeast of the TP and the cyclone over the west TP are the key circulation systems that dominate the influence of interdecadal Silk Road pattern on summer ITP rainfall.

The decadal variations of SPG SST and the interdecadal Silk Road pattern are well predicted in the ensemble mean of DCPP models, especially for the two-key upper troposphere circulation systems surrounding the TP (Fig. 4, E and F). Hence, we see a high prediction skill of the decadal variability of the upper troposphere winds over the TP (fig. S6), which leads to the skillful prediction of summer rainfall over the ITP. In addition, we find that the CMIP6 historical simulations can also partly reproduce the decadal variation of the summer ITP rainfall after 1980 (fig. S7), which is related to the natural forcing in driving the decadal variations of SPG SSTA after 1980 (39). This implies that the external radiative forcing could play a role in the predictability of summer ITP rainfall through modulating the SPG SST. In summary, the decadal predictability of the SSTA over the SPG regions in the North Atlantic and the associated interdecadal changes of the Silk Road pattern are covarying with the skillful decadal predictions of the summer ITP rainfall.

Real-time forecasts

How will the ITP rainfall change in the coming decade? Real-time predictions help policy-makers in adopting timely and effective mitigation and adaptation activities (18). We analyzed the 2- to 9-year mean forecast results of DCPP models with real-time forecast data available. In total, 60 ensemble prediction members are used (see Materials and Methods). We see a notable increase in rainfall over the large bodies of the TP in the near future (Fig. 5, A and B). The variance-adjusted forecast shows that the rainfall will increase in the central to the eastern TP but decrease in parts of the southwestern TP, with an amplitude of 0.10 to 0.30 mm/day. For the forecast of the summer rainfall over the ITP, the postprocessing procedure of variance adjustment can narrow down 41.8% uncertainties compared to the raw forecasts (see Materials and Methods), with the optimum estimate of 0.27 mm/day (0.11 to 0.41 mm/day; 5 to 95% range) during 2020–2027 relative to 1986–2005 (Fig. 5, C and D), indicating a 12.8% increase in summer mean precipitation.

Fig. 5 The 2020–2027 averaged summer TP rainfall anomalies in real-time forecasts.

(A and B) Spatial distributions of the 2020–2027 averaged summer precipitation anomalies (mm/day) from (A) decadal forecast from 2018 (averaged over raw forecast ensembles) and (B) decadal forecast from 2018 (averaged over variance-adjusted forecast ensembles). Regions with ACC < 0 are masked. (C and D) Probability distributions (%) of the ensembles for ITP summer precipitation anomalies from (C) raw forecast ensembles and (D) variance-adjusted forecast ensembles. The blue lines and numbers in (C and D) show the ensemble mean values, and the member sizes are also noted for each scenario. The precipitation anomalies are relative to the period of 1986–2005.


Numerous lakes in the interior ITP are affected by the decadal variability of summer ITP rainfall. For instance, marked lake shrinkages occurred on the ITP during the period of 1970–1990, but then, the lakes expanded significantly since the late 2000s, which play a vital role in regulating the water resources in the surrounding regions (26, 32). In addition, the hydroclimate responses of the TP lakes to rainfall have caused serial economical and environmental problems, such as desertification in the beach area and loss of grazing grassland, soil salinization and swamping (40), and transport infrastructure disruptions (41). Hence, a reliable climate prediction, especially several years ahead, is both crucial and urgent for the TP region. Our results will provide a useful reference to policy-makers for climate change adaptation and mitigation activities in these high-elevation regions that are sensitive to climate change.


Observational data

The following observational and reanalysis datasets are used: (i) The Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE; 1951–2015) (42), (ii) the GPCC monthly precipitation dataset (GPCC Precipitation 1.0 degree V2018 Full Reanalysis; 1891–2016) (43), (iii) the Global Precipitation Climatology Project (version 2.3; 1979–2020) (44), (iv) the Japanese 55-year Reanalysis Project (JRA-55, 1958–2017) (45), and (v) the Extended Reconstructed Sea Surface Temperature version-5 (ERSST v5, 1854–2017) (46).

Historical simulations and projections in CMIP6

We used monthly precipitation data from 19 CMIP6 models (table S1) for the projections of TP summer precipitation in Fig. 1B. Four realizations of projections under the Shared Socioeconomic Pathways (SSPs) SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 emission scenarios, respectively, and one realization of historical simulation are used for each model.

Decadal prediction experiments

The decadal prediction experiments used in this study are derived from the tier 1 set of experiments in the DCPP of CMIP6, with 138 ensemble members started from October, November, or December every year over the period of 1960–2016. To investigate the added value of initialization, we also use the historical simulations of the corresponding models from the CMIP6, which are independent of the initial conditions, consisting of 207 ensemble members over the period of 1950–2014. The multimodel ensemble mean is calculated by the unweighted mean of all ensemble members. In addition, 6 of the 10 climate prediction models used in our studies (MIROC6, EC-Earth3, NorCPM1, CMCC-CM2-SR5, HadGEM3-GC31-MM, and MPI-ESM1-2-HR) have real-time prediction results, with 60 ensemble members started from 2018, which also provide us with an opportunity to investigate the future changes in TP climate states. The analysis of atmospheric circulations is based on 9 of the 10 models excluding CMCC-CM2-SR5 because the hindcast of the monthly circulation data in CMCC-CM2-SR5 is unavailable. More details of the experiment used in this study are shown in table S2.

Prediction skill

We used the ACC and the MSSS to evaluate the deterministic prediction skill, that is, the ensemble mean prediction, which is an attempt to predict the most likely outcome and maximize the forecast skill (17).

The ACC and MSSS are defined asACC=i=1N(fif¯)(oio¯)i=1N(fif¯)2i=1N(oio¯)2(1)MSSS=1i=1N(fioi)2i=1N(oio¯)2(2)where N is the number of hindcast start dates and fi and oi are the ensemble mean prediction and the corresponding observation at each hindcast start date; overbar represents the average over all hindcast start dates. The ACC measures the prediction skill of the phase of variability, and the MSSS measures the prediction skill of the amplitude of variability.

For the decadal predictions, the predicted anomalies are calculated with respect to the predicted climatology, which is calculated by using hindcasts starting each year from 1960 to 2016. By using the anomalies, the mean state biases of the prediction for each start month and lead time are approximately removed, which is equivalent to the empirical bias correction. To be consistent with the predictions, the anomalies of the observed variables are calculated relative to the climatology of 1961–2016. The lead time of the prediction represents the time interval between the predictive target and the initial date. The 2- to 9-year averaged forecast of each start date is calculated by the rolling 8-year average of the forecast anomalies of 2 to 9 years. Similarly, the 4-year mean forecasts of 1 to 4 years, 2 to 5 years, 3 to 6 years, 4 to 7 years, 5 to 8 years, 6 to 9 years, and 7 to 10 years are calculated by the rolling 4-year average of the forecast anomalies of the corresponding forecast years. The 8-year (4-year) running averages of the observations are used to verify the rolling 8-year (4-year) averaged forecasts. Both observations and models are interpolated to 1.0° (longitude) × 1.0° (latitude) before comparison.

Contributions of initialized and uninitialized components to overall skill

A recently developed approach (19) is applied to investigate the relative contributions of the external radiative forcing and initialization to the overall decadal prediction skill of summer ITP rainfall. First, the time series of the observation (X), the decadal prediction (Y), and the historical simulation (U) are decomposed asX=χf+χ+x(3)YA={Yk}=ψf+ψ+{yk}(4)UA={Uk}=φf+{uk}(5)where k represents the ensemble members; χf, ψf, and φf represent the externally forced components; χ and ψ represent the predictable internally generated components; and x, yk, and uk represent the unpredictable noise components. The assumption is that the ψf and φf are common across the ensemble members, while the yk and uk are independent and identically distributed across their respective ensembles. When taking averages over large ensembles, which are denoted by the curly braces, the yk and uk are approximately equal to zero. The YA and UA are the ensemble mean of decadal prediction and historical simulation, respectively.

The uninitialized (initialized) components of the prediction is Yu (Yi), respectively, which areYu=αφf(6)Yi=(ψfαφf)+ψ(7)where α is the scalar projection of ψf on φf, which is defined asα=θrYAUAσYAσUAσφf2(8)where σYA, σUA, and σφf denote the SD of YA, UA, and φf, respectively; θ = 0 if the covariance between ψf and φf is negative, else θ = 1; and rYAUA is the correlation between the multimodel ensemble mean decadal predictions and the multimodel ensemble mean historical simulations.

The correlation skill can be decomposed asrXYA=rXYuσYuσYA+rXYiσYiσYA=ru+ri(9)where σYA, σYu, and σYi denote the variance of YA, Yu, and Yi, respectively; rXYA, rXYu, and rXYi are the correlation between the observation and ensemble mean decadal predictions, the external force, and the initialized components of the prediction, respectively; and ru and ri are the contributions of the initialized and uninitialized components to rXYA. The percentage contributions of the initialized and uninitialized components to rXYA are calculated by (rurXYA)×100% and (rirXYA)×100%, respectively.

Ratio of predictable component

The RPC (29) is defined asRPC=rXYA2rY(k)Yk2(10)where rXYA is the correlation between the multimodel ensemble mean prediction and the observation, rY(k)Yk is the correlation between the multimodel ensemble mean prediction and a single ensemble member, and the ensemble member is eliminated for the ensemble mean. The RPC is calculated by averaging over all possible combinations of ensemble members. For a perfect forecast system, the RPC should be equal to one. If RPC > 1, then the model can predict the real world better than itself (29, 47, 48), which is referred to as the signal-to-noise paradox (31, 47, 48).

Variance adjustment

A postprocessing procedure of variance adjustment (29, 30, 48, 49) is applied to make the amplitude of the predicted signal consistent with the observation. Unlike previous studies such as that of Eade et al. (29), where variance adjustment was based on the detrended time series, we keep the trend in the variance adjustment to realize the real-time forecast. Since the climate change signal can limit the effectiveness of adjustment if the climate change signal expresses as a long-term trend (see text S1), the components of long-term linear trend are considered in the method of variance adjustment used in this study.

First, the time series of the observation and the decadal prediction are decomposed asX=XXtrend(11)YA=YAYAtrend(12)Yk=YkYAtrend(13)where X′, YA, and Yk represent the detrended value of observation, ensemble mean prediction, and ensemble member prediction and the long-term linear trend time series of the observation and the original hindcast ensemble mean are Xtrend and YAtrend, respectively.

The variance adjustment to the predictable components of the raw prediction is as followsYA˜=YA·σXR(X,YA)σYA+YAtrend·σXtrendσYAtrend(14)where YA is the adjusted hindcast ensemble mean value, σX (σYA) is the SD of X′ and YA, R(X,YA) is the correlation coefficient between X′ and YA, and σXtrendYAtrend) is the SD of the Xtrend and YAtrend.

The variance adjustment to the ensemble members’ variance about the adjusted ensemble mean is as followsYk˜=Yk·σX[1R2(X,YA)]σnoi+YA˜(15)where Yk is the adjusted ensemble member value and σnoi is the SD of the noise, which is written as followsσnoi=variance(YkYA)=variance(YkYA)(16)

For the real-time predictions, the uncertainty is defined as the range between the minimum and maximum values.

Eddy geopotential height

To verify the mid-latitude wave-like pattern more clearly, the eddy geopotential height anomalies (50) are shown in Fig. 4 and fig. S5, which are defined as the deviation of the height anomalies from the zonal mean.

Wave activity flux

We use the WAF (51) to measure the wave propagation. The WAF represents the horizontal propagation of the quasi-stationary Rossby waves. Its horizontal components in pressure coordinates are as followsWAF=p2u¯{u¯(ψx2ψψxx)+v¯(ψxψyψψxy)u¯(ψxψyψψxy)+v¯(ψy2ψψyy)}(17)

Here, overbars and primes denote mean states and low pass–filtered anomalies, respectively. The quantity u¯ denotes the magnitude of the climatological wind. u and v represent the zonal and meridional wind, respectively. ψ is the stream function, and p is the normalized pressure (pressure/1000 hPa). The subscripts x and y represent the zonal and meridional gradients, respectively.

Maximum covariance analysis

The Silk Road pattern can be defined as the first empirical orthogonal function mode of the summer mean 200-hPa meridional wind over the Eurasian continent (52). To extract the coordinated interdecadal variation of SSTA over the SPG region and the atmospheric disturbances over the downstream regions, we use the MCA (53, 54) to the SSTA over the SPG region (50° to 65°N, 60° to 10°W; Fig. 4D, blue box) and the summer 200-hPa meridional wind anomalies over the mid-latitude regions of Eurasia (20° to 70°N, 50°W to 100°E; Fig. 4C, blue box). For the first MCA mode in the observation, the value of squared covariance fraction is 86%, and the correlation coefficient between two expansion coefficients of the first MCA is 0.89, indicating that the SSTAs over SPG region and Silk Road pattern are tightly coupled on the interdecadal time scale.

Statistical analysis

The significance levels of the skill (ACC or MSSS) and RPC are tested by a nonparametric bootstrap approach (30, 33, 55), as shown in Figs. 2 and 4 (E and F). For uncertainties of the values, a finite ensemble size (E) and a finite number of validation points (N) are taken into account. The detailed steps are as follows:

1) Randomly sample with replacement N validation cases. To take autocorrelation into account, every five consecutive years of the observations and predictions are composed as a block to conduct resampling and replacement.

2) For each case, randomly sample with replacement E ensemble members.

3) Compute the statistic for the ensemble mean (e.g., ACC, MSSS, or RPC).

4) Repeat (1) to (3) 1000 times to create a probability distribution.

5) Obtain the significance level based on a two-tailed test of the hypothesis that skill (ACC or MSSS) is zero or RPC is one.

For skill as a function of ensemble size, we randomly sample without replacement to obtain the required number of ensemble members and compute the average for all possible combinations. Uncertainties in Fig. 3D are computed on the basis of single members because the samples are not independent for larger ensembles. To test the significance of regression analysis, Student’s t test was used in Fig. 4 and fig.S5.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge the World Climate Research Programme’s Working Group on Coupled Modeling, which is responsible for the CMIP6 climate modeling groups and the decadal prediction project (listed in tables S1 and S2) for producing and making available their model output. We thank W. Zhang for helpful comments and discussions. Funding: The study was supported by the Chinese Academy of Sciences under grant no. XDA20060102, the National Natural Science Foundation of China under grant no. 41988101, and the International Partnership Program of Chinese Academy of Sciences under grant no. 134111KYSB20160031. Author contributions: T.Z. designed the research, provided comments, and revised the manuscript. S.H. performed the analysis and drafted the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The historical climate simulations from CMIP6 listed in table S1 and the decadal prediction experiments from CMIP6 DCPP listed in table S2 are available at The observational precipitation and SST datasets are available on the Climate Data Guide website ( The Japanese 55-year reanalysis dataset is available at The data in this study were analyzed, and the figures were created with NCAR Command Language (NCL). All relevant codes used in this work are available from

Stay Connected to Science Advances

Navigate This Article