Colorado River water supply is predictable on multi-year timescales owing to long-term ocean memory

Skillful multi-year climate forecasts provide crucial information for decision-makers and resource managers to mitigate water scarcity, yet such forecasts remain challenging due to unpredictable weather noise and the lack of dynamical model capability. Here we demonstrate that the annual water supply of the Colorado River is predictable up to several years in advance by a drift-free decadal climate prediction system using a fully coupled climate model. Observational analyses and model experiments show that prolonged shortages of water supply in the Colorado River are significantly linked to sea surface temperature precursors including tropical Pacific cooling, North Pacific warming, and southern tropical Atlantic warming. In the Colorado River basin, the water deficits can reduce crop yield and increase wildfire potential. Thus, a multi-year prediction of severe water shortages in the Colorado River basin could be useful as an early indicator of subsequent agricultural loss and wildfire risk. Sea surface temperature anomalies in the Atlantic and Pacific Oceans can help predict water shortages in the Colorado River basin, according to analyses of decadal climate predictions and observations.

T he Colorado River is the most important water resource in the semi-arid western United States (U.S.). Demand for its water has increased continuously for the past 60 years but its water supply is facing an unsustainable future due to declining precipitation and prolonged drought events [1][2][3] . Prolonged water shortages in the Colorado River can cause serious damages to a wide range of sectors, including agriculture, forestry, energy, food security, drinking water, and tourism 4,5 . High-impact drought events in the Colorado River basin during 2000-2016 severely stressed regional water supply and strained reservoir operations, recreation, and ecological services 1,3 . To help managers and policy-makers cope with drought-induced water shortages coupled with ever-increasing demand, it is crucial to develop multiyear drought predictions for safeguarding and maintaining industrial and societal wellbeing.
Presently, U.S. operational drought forecasts primarily focus on day-to-month outlooks, such as the short-term drought severity indicator and the monthly drought outlook generated by the Climate Prediction Center (CPC) and the National Drought Mitigation Center of the National Oceanic and Atmospheric Administration. These forecasts are severely limited by shortterm weather phenomena 6 , since unpredictable atmospheric noise makes interannual-to-decadal climate prediction challenging 7,8 . However, previous research has identified a prominent "drought cycle" in the Intermountain West and the southwest U.S., which was attributed to climate variability associated with long-term ocean memory [9][10][11][12][13][14] . Besides ocean memory, land systems (i.e., soils, groundwater, streamflow, vegetation, and perennial snowpack) filter out the high-frequency precipitation fluctuation and integrate atmospheric signals over space and time [15][16][17][18][19] . These results lead to a hypothesis proposed herein that skillful multi-year predictions of the Colorado River water supply are possible by utilizing long-term ocean memory, its associated atmospheric teleconnections, and the natural filtering effect in the land system altogether. A similar concept has been applied to develop seasonal drought forecasts based on El Niño Southern Oscillation (ENSO) predictions 20 and a statistical model for multi-year water supply predictions 21,22 , yet its application to the water supply forecast beyond seasonal timescales remains unknown due to the signal-to-noise paradox 7 .
Here we demonstrate that interannual-to-decadal variability of the Colorado River water supply is predictable for several years in advance using a decadal climate prediction approach. Our assessment is based on three experiments using the fully coupled climate model Community Earth System Model (CESM; see "Methods"): a historical and future emission scenario simulation (the externally forced run), an ocean data assimilation run (referred to as the ASSIM run), and a multi-year initialized prediction experiment (the hindcast run). In the ASSIM run, we assimilated the observation-based 3-dimensional ocean temperature and salinity anomalies 23 into the ocean component of CESM, with prescribed natural and anthropogenic radiative forcings. As a result, the model-simulated atmosphere-land variability in the ASSIM run is indicative of the response to ocean variability and external radiative forcings. We conducted the ocean assimilation experiment for the period 1960-2015 with 10 ensemble members. To minimize the bias known as artificial model drift during the simulation, an effective bias-adjustment method was adopted and applied in the assimilation method 24 . Predictability of Colorado River water supply and its sources are determined by the externally forced run and the hindcast run.

Results
Reconstruction of Colorado River water supply. To detect water shortages of the Colorado River, we used a historical record of water supply in the Colorado River basin provided by the Bureau of Reclamation 1 (blue line in Fig. 1a). This Colorado River water supply dataset was designed to develop adaptation and mitigation strategies for water resource agencies and stakeholders throughout the Colorado River basin. Its interannual-to-decadal variability mostly reflects the natural flow at Lees Ferry, Arizona, and is assumed to be free of human water usage (e.g., irrigation) 25 . This water supply record illustrates severe shortages in the years 1963,1977,1981,1990,2002,2012, and 2013 (yellow shade in Fig. 1a). These severe shortages correspond with agricultural losses and high fire activity in the Colorado River basin, as described later.
According to previous studies, the major fraction of Colorado River streamflow variability arises from groundwater variability through precipitation input as well as land processes of infiltration, subsurface storage and transmission, and convergence toward the channels 17,26,27 . Through these regional processes, Colorado River streamflow is considered to have a large spatial footprint of hydroclimate conditions in the Intermountain West 17 . To reveal regional climatic control of the Colorado River water supply, we made correlation maps of precipitation and total soil water (i.e., the sum of soil water in all soil layers) as observation-based and model-simulated anomalies with the observed Colorado River water supply ( Fig. 1b-e). While the correlation between the Colorado River water supply and annual mean precipitation variability is high, we find more significant correlations of water supply with soil water anomalies across the western U.S. (Fig. 1b, d). The ASSIM run also demonstrates that the Colorado River water supply positively correlates with soil water anomalies around Utah, Colorado, and New Mexico (Fig. 1e), albeit with somewhat distorted spatial patterns of precipitation and soil water anomalies compared to the observation. By taking the regional average over these highly correlated soil water anomalies (black boxes in Fig. 1b, e), we reconstruct the Colorado River water supply using the National Centers for Environmental Prediction CPC reanalysis and the ASSIM run (black lines in Fig. 1a). This reconstructed water supply shows a close agreement with the observed temporal variation of the Colorado River water supply. This result supports our hypothesis that the Colorado River water supply is closely linked to soil water variability in the Intermountain West.
Atmospheric weather disturbances still generate unpredictable high-frequency noise, inhibiting precipitation predictability, yet such high-frequency components of precipitation variability are mostly filtered out through the land hydrological processes in the Intermountain West 11,15 . Because of this land filtering effect, annual soil water anomalies reflect precipitation anomalies averaged for several years ( Supplementary Fig. 1) that link with long-term ocean memories. Consequently, we reconstructed the Colorado River water supply from the areal average of soil water anomalies associated with interannual-to-decadal ocean variability in the ASSIM run. The 10-member ensemble mean of ASSIM runs demonstrates skill in capturing many major historical shortages of the Colorado River water supply (Fig. 1a), even though it did not include any atmospheric or land observations. The correlation coefficient between the observed and the model reconstructed Colorado River water supply was higher when we applied a 3-year running mean filter to highlight the low-frequency climate variability (correlation coefficients are 0.43 for annual and 0.60 for 3-year means, respectively). Our analysis suggests that the Colorado River water supply is potentially predictable through "perfect knowledge" of ocean conditions. Predictability of Colorado River water supply. Next, we examined the dynamical retrospective forecasts for the period 1960-2015, which was initialized annually on January 1st based on the ASSIM run (referred to as the hindcast run; see Methods). This hindcast experiment consists of 10-year-long predictions with 10 ensemble members for each initialized run. The ensemble mean of the hindcast runs for the reconstructed Colorado River water supply and the corresponding ensemble spread aligned well with the CPC reanalysis and the ASSIM run for 1-year and even 2-year lead times (Fig. 2a, c). The predictive skills measured by the anomaly correlation coefficient (ACC) and the root-meansquare error (RMSE) in our hindcast run outperform those in the externally forced run and persistent forecast (see "Methods"), up to the lead time of 48 months (Fig. 2b, d). These higher skills in the hindcast run indicate that ocean initialization is crucial for skillful forecasts of the Colorado River water supply on interannual-to-decadal timescales whereas atmospheric initialization can improve land hydrological predictability on seasonal Observed and reconstructed Colorado River water supply. a Annual mean time series of Colorado River water supply (blue) and its reconstructions by CPC soil moisture reanalysis (black broken) and the ensemble mean of 10-member ASSIM runs (black solid). The reconstructions of the Colorado River water supply are obtained by area averaged soil moisture anomalies in CPC reanalysis (33°N-43°N, 119°W-103°W) and ASSIM runs (black boxes in d, e). Years with prominent water supply shortages are highlighted with yellow shading (1963, 1977, 1981, 1990, 2002, 2012, and 2013). Correlation maps of b, c precipitation and e, f soil moisture anomalies in the observation-based estimate (left) and the ASSIM run (right) against the Colorado River water supply for 1960-2015. Red line and yellow mark in b, d indicate geographical locations of Colorado River and Lees Ferry, respectively. Correlation coefficients of 0.28, 0.33, 0.45, and 0.53 correspond to the statistical significance at 90%, 95%, 99%, and 99.9% levels with 36 degrees of freedom on the basis of twosided Student's t-tests and an equivalent sampling size of the Colorado River water supply. Anomalies are defined as deviations from climatological means and linear trends are removed in each grid. timescales 28 . As a result, the hindcast run demonstrates skill in predicting the Colorado River water supply for 2 years in advance (Fig. 2b, d). Furthermore, we find higher predictive skill by measuring the hindcast run against the model reconstruction of the ASSIM run (i.e., red dotted lines in Fig. 2b, d) compared to either the observation-based products of CPC reanalysis or the observed Colorado River water supply. Because of deficiencies in current climate models simulating observed land hydrology, these results suggest that enhanced predictive skill could be achieved by improving model performance in simulating soil water variability. The dynamical prediction presented here supports earlier statistical predictions of Colorado River water supply years ahead 21,22,29,30 .
Our prediction system suggests that the changes in large-scale atmospheric circulations link with water shortages of the Colorado River. During extreme water shortages, observations show drier soil water conditions around the Intermountain West under atmospheric high-pressure anomalies with a North Pacific connection (Fig. 3g). The drier soil patterns appear one year before the water shortages and are limited to the eastern side of the Intermountain West (Fig. 3d). We also see higher and lower atmospheric pressure anomalies in the North Pacific 1 and 2 years before the shortages (Fig. 3a, d), implying that the water shortages are induced by an atmospheric teleconnection, such as the Pacific North American (PNA) pattern. However, the statistical significance of these anomalies is obscured by highfrequency atmospheric noise. These observed PNA-like atmospheric patterns and drier soil conditions in the Intermountain West are well captured in the ASSIM run 1 and 2 years before water shortages (center panels in Fig. 3). The 1-year prediction in the hindcast run also captures these PNA-like atmospheric teleconnection patterns and dry soil water conditions (right panels in Fig. 3). These results strongly support the hypothesis that the annual Colorado River water supply is predictable by utilizing ocean memory and the land filtering effect.
It is worth noting discrepancies between the observations and the model simulations. As shown in Fig. 3g, e, the ASSIM run exhibits the timing of water shortages a year too early. In addition Fig. 2 Predictability of reconstructed Colorado River water supply. Time series of reconstructed Colorado River water supply and the hindcast runs for a 1-12 and c 13-24 months lead time. Solid and broken black lines indicate reconstructions by the ASSIM run and the CPC reanalysis, respectively. Red circles and bars denote ensemble mean and spread (±1 standard deviation) of 10-member hindcast run initialized on 1st January, every year. Blue lines and shading denote ensemble mean and spread (±1 standard deviation) of the 10-member externally forced run. Anomalies are defined as deviations from climatological means and linear trends are removed in each grid. Time series are normalized by one standard deviation based on the CPC reanalysis and the ocean data assimilation run. b, d Predictive skills of reconstructed Colorado River water supply in hindcast run (red solid), externally forced run (blue), and persistence (black) measured by anomaly correlation coefficient (ACC) and root-mean-square-errors (RMSEs) against the CPC reanalysis (i.e., black dashed lines in a, c) based on 50 initialized cases (1961-2010). A 12-month running mean filter is applied to linearly detrended anomalies. Red stars in b, d denote predictive skills in the hindcast run against the observed Colorado River water supply (i.e., blue line in Fig. 1a). Red dotted lines in b, d correspond to the predictive skills in the hindcast run against the ocean data assimilation run (i.e., black solid lines in a, c), which indicates the upper limit of predictability in the current model (i.e., potential predictability). In the persistent predictions, the observed anomalies in initial years (i.e., time average of anomalies for 0-11 months before the start time of prediction) are maintained during the prediction lead time. Black broken lines in b, d denote a correlation coefficient of 0.235 with the statistical significance at 95% level with 50 degrees of freedom using one-side Student's t-test and one standard deviation. Fig. 3 Atmospheric teleconnections associated with water shortages of the Colorado River. Composite maps of annual mean soil water (shaded) and 500 hPa geopotential height anomalies (contours) in observations (left), the ASSIM run (center), and the 1-12 months hindcast run (right) at (a-c) 2 years before, d-f 1 year before, and g-i during waters shortage in the Colorado River basin (i.e., yellow shading in Fig. 1). Anomalies are defined as deviations from climatological means and linear trends are removed in each grid. Contour intervals are ±3, ±6, ±9, ±12, ±15, ±18, and ±21 hPa. Zero contours are omitted and negative values are dashed. Light and dark pinks (greens) indicate the statistical significance at the 90% and 95 % levels with two-side Student's t-test for soil water (geopotential height) anomalies. COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-020-00027-0 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | (2020) 1:26 | https://doi.org/10.1038/s43247-020-00027-0 | www.nature.com/commsenv to this earlier timing, the centers of action in the PNA-like atmospheric teleconnection are slightly shifted in the ASSIM run compared to the observations (Fig. 3). It is also unclear what the model capability is in simulating evapotranspiration that plays an important role in the land filtering effect. Because of these model deficiencies, we selected a different area for reconstructing the Colorado River water supply in the ASSIM run compared to the CPC reanalysis (boxes in Fig. 1d, e). These model deficiencies may reduce skill in predicting observed Colorado River water supply in our hindcast run compared to the skill against the ASSIM run (start marks and dotted lines in Fig. 2b, d). Further enhancement of predictive skills could be achieved by improving model-simulated land hydrological processes and using finer model spatial resolution.
Ocean precursors for water shortage. Given the demonstrated ocean influences on multi-year predictability in the Colorado River water supply, in which drought years in the ASSIM run are linked to interannual-to-decadal ocean variability, we extract 7 predictable years when the Colorado River water supply is more than −0.5 standard deviations from the mean in the observed dataset and the reconstruction of ASSIM run ( Supplementary  Fig. 2): 1963,1976,2000,2001,2002,2012, and 2013. The composite maps associated with these years exhibit the PNA-like atmospheric teleconnection pattern 1-3 years before the mature phase (Fig. 4). Significant changes in SST are also noticeable, such as a La Niña-like SST pattern in the tropical Pacific 1-2 years before, a warmer SST in the Kuroshio-Oyashio extension region 1-3 years before, and a warmer SST in the southern tropical Atlantic 2-3 years before the mature phase in both the observations and the ASSIM run (Fig. 4). We consider these SST anomalies precursors for severe water shortages of the Colorado River.
To evaluate ocean impacts on water shortage, we develop a statistical prediction based on the three SST precursors for the Colorado River water supply (Fig. 5). To capture the lowerfrequency component, we apply a 2-year running mean to the area-averaged SST anomalies. Lead-lag relationships reveal that the observed Colorado River water supply correlates with SST precursors in the tropical Pacific at 1-2 years lag, the North Pacific at 2-3 years lag, and the southern tropical Atlantic at 3-4 years lag (Fig. 5b). A multiple regression analysis using these inputs indicates that this statistical prediction model can explaiñ 38% of the variance in Colorado River water supply (correlation coefficient R = 0.62; Fig. 5a). Specifically, prominent negative values of this statistical model (less than −1.0 standard deviation) always accompany drier anomalies of the Colorado River water supply. In other words, by monitoring these SST precursors, we can detect emergence of water shortage in the Colorado River basin. The contributions from each precursor correspond to 7% from the tropical Pacific, 34% from the North Pacific, and 59% from the southern tropical Atlantic (see Methods). According to this estimate and the different lead times of SST precursors, it may be possible to provide an early indicator of water shortage at least two years before a drought event by utilizing the North Pacific and the southern tropical Atlantic SST precursors. Further improving predictive skills could include developing a more sophisticated statistical method using three ocean precursors as well as combining the dynamical predictions of these ocean precursors with the statistically-based prediction model.
Drought impacts on agriculture and wildfire. The oceaninduced large-scale climatic anomalies modulating the Colorado River water supply also affect crops and natural ecosystems. Even though advances in agricultural technology and improved management practices have led to a remarkable increase in annual crop yields (Fig. 6a, since the early 1960s 31 ), crop production remains vulnerable from prolonged drought events. During severe water shortages (yellow shade in Figs. 1 and 6), we see large reductions in annual crop yields for most major crops (alfalfa, wheat, corn silage, and barley; Fig. 6a) in the Intermountain West (encompassing Colorado, Nevada, Utah, Wyoming, and Idaho). Specifically, small grains (wheat and barley) are more vulnerable to drought stress than alfalfa and corn silage (Fig. 6a and Supplementary Table 1). The average across linearly detrended and normalized crop yields further highlight the effect of these prolonged droughts ( Fig. 6b and Supplementary Fig. 3). Interannual-to-decadal variability in the averaged crop yield significantly correlates with the Colorado River water supply (R = 0.59). Consistent with this positive correlation, we can see a shift in the probability distribution toward Mexico. The observed Colorado River water supply and its reconstruction in the assimilation run in Fig. 1a are replotted in b, c as well (blue and black lines). Years with prominent water supply shortages are highlighted by yellow shading (1963, 1977, 1981, 1990, 2002, 2012, and 2013). Annual crop yields are represented by a percentage relative to the values in 1960. The crop yield in b was linearly detrended in each state, each crop, averaged across five states, normalized by one standard deviation, and then we took the ensemble mean of all five crops (see Supplementary Fig. 3). Annual mean fire size is calculated by fire burned area divided by the number of large fire events per year.
reductions in crop yield during shortages of the Colorado River water supply (Fig. 7a).
Severe water shortages also influence the occurrence of large wildfires 5,32 throughout the West (Fig. 6c; including Colorado, Nevada, Utah, Wyoming, Idaho, Oregon, California, Arizona, and New Mexico). Whereas interannual-to-decadal wildfire variability demonstrates complex behavior by interacting with soil water and vegetation 15 , drought is one of the main factors  33 . We found that interannual-to-decadal variability of mean fire size positively correlates with the Colorado River water supply (R = 0.47, Fig. 6c). Large fire events occurred during the severe water shortages in 1990, 2002, 2012, and 2013, which is also highlighted by an increase in the probability distribution of large fires (Fig. 7b). Since water shortages of the Colorado River are tightly linked with drier soil conditions across the basin (Figs. 1d and  7c), our multi-year predictions of Colorado River water supply can apply to risk assessments of crop yield reduction and large fire occurrence.

Discussion
Using the decadal climate prediction system demonstrated here, we found that the Colorado River water supply is potentially predictable by utilizing long-term ocean memory and the land filtering effect. However, sources of long-term ocean memory remain unclear. Many previous studies have pointed to the linkage of drought conditions in the U.S. with interannual-todecadal sea surface and climate variability in the tropical and North Pacific, such as ENSO 34 , Pacific Decadal Oscillation 35 , Interdecadal Pacific Oscillation 36 , and Pacific Quasi-decadal Oscillation 11,14 . The SST anomalies associated with these phenomena serve as an atmospheric heat source and so, induce a PNA-like atmospheric response through the atmospheric Rossby wave propagation. This PNA-like atmospheric pattern can persist as a response to lower-frequency ocean forcing while modulating moisture transport from the ocean to the Colorado River basin. Recent studies have also introduced the concept of inter-basin climate interactions 37,38 , in which the Atlantic Ocean can modulate Pacific climate variability [39][40][41] . The demonstrated predictions of the Colorado River water supply are in good agreement with this concept. Whereas some previous studies have pointed out the relationship between U.S. drought conditions and the North Atlantic SST anomalies associated with the Atlantic Multidecadal Oscillation 22 , our analysis detects a more important role of the southern tropical Atlantic.
The prediction of the Colorado River water supply obtained here has implications not only for drought management but also for proactive management practices that might be applied to mitigate crop yield reductions as well as wildfire risk (Figs. 1 and  6). Our results indicated that multi-year drought conditions in the Intermountain West are remotely linked to three ocean precursors, including the La Niña-like SST cooling in the tropical Pacific 1-2 years before, the warmer SST in the Kuroshio-Oyashio extension region 2-3 years before, and the warmer SST in the southern tropical Atlantic 3-4 years before the mature stage of Colorado River water shortages. By monitoring these ocean precursors, we can develop a trustworthy outlook for heightened drought threat in the upcoming year. The lead time could help water managers and stakeholders establish mitigation strategies and proactive management plans, such as deliveries, allocations, conservation, and efficient usages of operational water supply on an annual basis. The long-range drought outlook is also applicable to inform farmers, ranchers, forest and range managers, energy supply potential, and fire managers about the potential for upcoming water scarcity and subsequent risks. This approach of using optimized model-simulated soil water variability to predict multi-year water supplies of the Colorado River basin could also be used for other major river basins in the western U.S. and may have other implications for energy supply potential in large hydropower plants or socioeconomic risks.

Methods
Observation data. Annual Colorado River water supply was obtained from a comprehensive Colorado River basin Water Supply and Demand Study by the Bureau of Reclamation 1 . Annual crop yield for five dominant crops (in terms of land area) in the Mountain West (alfalfa, spring wheat, winter wheat, corn silage, and barley) for 1960-2015 were obtained from agricultural census and survey data collected by United States Department of Agriculture-National Agricultural Statistics Service 42 . State-level mean yields for Colorado, Nevada, Utah, Wyoming, and Idaho were averaged together by year. Annual mean wildfire size for 1984-2014 was calculated by dividing wildfire area burned by wildfire occurrences in Colorado, Nevada, Utah, Wyoming, Idaho, Oregon, California, Arizona, and New Mexico. These observational wildfire datasets were derived from the Monitoring Trends in Burn Severity database 43 , which includes fires >405 ha in the western US. Climate datasets for 1960-2018, such as precipitation, soil water, sea surface temperature, and geopotential height at 500 hPa, were derived from gridded datasets of the Global Precipitation Climatology Centre (GPCC) combined full version 7 and version 4 monitoring data products 44 , the National Centers for Environmental Prediction Climate Prediction Center (CPC) reanalysis product 45 , Characteristics of Global Sea Surface Temperature Analysis Data version 2 46 , and Japanese 55-year reanalysis datasets 47 , respectively.
Model experiment. The fully coupled climate model used here is a lowerresolution version of the Community Earth System Model (CESM) version 1.0.6 48 . The atmospheric and land models correspond to a horizontal resolution of T31 spectral grid (approximately 3.75°horizontal resolution) with 26 hybrid sigma/ pressure coordinate levels in the atmosphere, 10 soil layers, and an aquifer water layer on land. The land component is the Community Land Model version 4 (CLM4), which includes a carbon-nitrogen biogeochemical cycle, a simple groundwater model, and a wildfire scheme 49 . The ocean and sea-ice models use a curvature grid with a displaced North Pole (~3°horizontal grid but 1°latitude grid near the equator) with 60 vertical levels in the ocean. This computationally efficient model resolution enables the running of a large number of model integrations for the decadal climate prediction experiment. No flux correction is applied in exchanging heat, water, and momentum fluxes between the atmosphere and the ocean. Details of the basic model performance in this configuration can be found in previous studies 15,24,48 .
Using CESM, we conducted three experiments based on the decadal climate prediction protocol: an externally forced run, an ocean data assimilation run, and a series of multi-year hindcast runs. These experiments are based on a drift-free decadal climate prediction system 24 . All experiments were conducted with a 10member ensemble and use the same radiative forcings (i.e., greenhouse gas and aerosol concentrations, solar cycle variations, and major volcanic eruptions) as the historical record before 2005 and the IPCC RCP4.5 future emissions scenario after 2005 50 . In the externally forced run, the model was prescribed by natural and anthropogenic radiative forcings for the period 1850-2030. Initial conditions for the externally forced run were obtained from 10 random years of the pre-industrial control simulation using a constant external forcing for the year 1850. From this externally forced run, we obtained 10-member initial conditions on 1st January, 1958 and assimilated the observation-based 3-dimensional ocean temperature and salinity anomalies into the ocean component of CESM for the period 1958-2015 (the ASSIM run). Based on the initial conditions obtained from the ASSIM run, we conducted a series of 10-year-long hindcast/forecast experiments from 1960 to 2015 every year.
To minimize artificial model climate drift during prediction, we assimilated the observed internal variability (Y int ) while maintaining the model simulated states in the climatology (X clm ) and the externally forced signal (X ext ) (i.e., the bias-adjusted observations is equal to Y int + X ext + X clm ), where Y is the observation, X is the ensemble mean of the externally forced run, and subscripts int, ext, and clm are internal variability, externally forced signal, and climatology, respectively. As a result of applying the balanced ocean conditions to the model-simulated climatological and externally forced signals, our hindcast run shows only little drift and no initialization shocks. This new drift-free prediction system has exhibited strong multi-year predictive skill for decadal climate variations of the Atlantic meridional overturning circulation, North Pacific decadal variability, and droughtfire variability in southwestern North America 9,24 .
Predictive skills. Predictive skills in the hindcast run, the externally forced run, and the persistent forecast were measured based on 50 initial cases (i.e., every year from 1961 to 2010). The hindcast run is 10 years long and uses 10-member ensemble predictions for every initial case. To assess the impact of ocean initialization, we converted the 10-member externally forced run into the same format as the hindcast run. In these runs, anomalies were calculated by removing the Fig. 7 Impact of Colorado River water shortages on agriculture and wildfire. Probability distributions of a the normalized and linearly detrended annual crop yields in the Mountain West (barley, corn silage, spring wheat, winter wheat, and alfalfa in Colorado, Nevada, Utah, Wyoming, and Idaho), b mean fire size in the West (Colorado, Nevada, Utah, Wyoming, Idaho, Oregon, California, Arizona, and New Mexico), and c the linearly detrended soil water anomalies in the Colorado River basin during all years (black) and the years of Colorado River water supply shortages (red). Maximum sample sizes for all years correspond to 1400 (56 years × 5 crops × 5 states) for crop yield, 279 (31 years × 9 states) for fire size, and (56 years × 639 grid points) for soil water although actual sample sizes for crop yield and fire size (n) are smaller than these maximum sizes due to missing values. Seven years of Colorado River water supply shortage are shown in yellow shading in Fig. 1. The horizontal axis of mean fire size b corresponds to the log10 scale. climatology mean and linear trends for 50 years as a function of lead time. Their anomalies in Fig. 4 were normalized based on the standard deviation of the ASSIM run and were applied by a 12-month running mean filter. We labeled the lead time based on the last month of the averaged period (i.e., the lead time labeled "12 months" corresponds to the average of 1-12 months lead time), which means that we do not have values for lead times labeled 1-11 in the hindcast run. In the persistent forecasts, initial anomalies obtained from the CPC reanalysis persist during the forecast period.
A statistical model based on the SST precursors. We developed the statistical model using multiple regression of three SST precursors: SST anomalies averaged over the tropical Pacific (15°S-15°N, 150°W-120°W), the North Pacific (25°S-35°N, 130°E-180°), and the southern tropical Atlantic (20°S-5°N, 30°W-20°E). Anomalies are defined as deviations from the climatological mean and the linear trends are removed at each grid point. These area averaged SST anomalies are normalized by one standard deviation in each region (Supplementary Fig. 4). Based on the leadlag correlation, we selected lead times of 1-2 years lag in the tropical Pacific (TP), 2-3 years lag in the North Pacific (NP), and 3-4 years lag in the southern tropical Atlantic (TA). The area in the tropical Pacific SST precursor is optimized based on SST anomaly pattern at a 1-year lag (Fig. 4e, f). Correlation coefficients between three SST precursors are summarized in Supplementary Table 2. A multiple regression analysis for these precursors with the observed Colorado River water supply provides the prediction (y) as follow: y = 0.074 × TP−0.337 × NP −0.589 × TA.

Data availability
The Colorado River water supply data were provided by Dr. Lawrence Hipps at Utah State University (lawrence.hipps@usu.edu) and James Prairie at the Bureau of Reclamation (jprairie@usbr.gov). Observed sea surface temperature dataset and reanalysis product of geopotential height were provided by Japan Meteorological Agency at http://ds.data.jma.go.jp/tcc/tcc/products/elnino/cobesst/cobe-sst.html and https://jra. kishou.go.jp/JRA-55/index_en.html. CPC soil water data were provided by the NOAA/ OAR/ESRL PSL, Boulder, Colorado, USA, at https://www.psl.noaa.gov/data/gridded/ data.cpcsoil.html. Other datasets generated in this study are available at https://climate. usu.edu/people/yoshi/data/2020-CoRiv_data/data.html. All datasets are the NetCDF format. The datasets generated in this study and figures are also available on request from the corresponding author.

Code availability
We used the NCAR Command Language version 6.5.0 51 to analyze observed and modelsimulated data and draw figures. The source code of the fully coupled climate model CESM 1.0.6 is distributed by the National Center for Atmospheric Research (http://www. cesm.ucar.edu/models/cesm1.0/). Ocean data assimilation and visualization codes are available on request from the corresponding author.