Characterizing the multisectoral impacts of future global hydrologic variability

There is significant uncertainty in how global water supply will evolve in the future, due to uncertain climate, socioeconomic, and land use change drivers and variability of hydrologic processes. It is critical to characterize the potential impacts of uncertainty in future water supply given its importance for food and energy production. In this work, we introduce a framework that integrates stochastic hydrology and human-environmental systems to characterize uncertainty in future water supply and its multisector impacts. We develop a global stochastic watershed model and demonstrate that this model can generate a large ensemble of realizations of basin-scale runoff with global coverage that preserves the mean, variance, and spatial correlation of a historical benchmark. We couple this model with a well-known human-environmental systems model to explore the impacts of runoff variability on the water and agricultural sectors across spatial scales. We find that the impacts of future hydrologic variability vary across sectors and regions. Impacts are felt most strongly in the water and agricultural sectors for basins that are expected to have unsustainable water use in the future, such as the Indus River basin. For this basin, we find that the variability in future irrigation water withdrawals and irrigated cropland increase over time due to uncertainty in renewable water supply. We also use the Indus basin to show how our stochastic ensemble can be leveraged to explore the global multisector consequences of local extreme runoff conditions. This work introduces a novel technique to explore the propagation of future hydrologic variability across human and natural systems and spatial scales.


Introduction
Water is critical for sustaining life and is central to nearly every human endeavor, including food and energy production [1].Future water availability across regions is important to characterize yet deeply uncertain, given the complex nature of the global hydrologic cycle and the dependence of water supplies and demands on highly uncertain climate and socioeconomic conditions [2][3][4].While previous studies have explored the impacts of uncertain future climate and socioeconomic trajectories on water availability and its impacts on the agricultural and energy sectors (e.g.[5][6][7]), the uncertainty inherent in the hydrologic system and its propagation across sectors has not been systematically explored at a global scale.The primary objective of this work is to explore how variability in projections of future surface water supply impact regional-toglobal food and water resource systems by (1) developing a novel stochastic model of future runoff with global coverage and (2) propagating the results of this stochastic model through a human-environmental system model with representations of energy, water, land, climate, and the economy.
To explore how variability in projections of future water supply impact regional-to-global food and water resource systems, we conduct our analysis using the global change analysis model (GCAM) [8], a multisector dynamics (MSD) model [9] that represents the interactions of the climate, energy, land use, and water systems with global coverage.MSD models have been developed across spatial scales (e.g.[9][10][11]) and commonly integrate multiple sectorspecific models, including climate, demographic, and hydrologic models.In the past, global multisectoral studies have often relied on global hydrologic models (GHMs) to explore plausible hydrologic futures at the global scale.GHMs are used to simulate the complex dynamics of the global hydrologic cycle and are therefore useful for quantifying global water availability and capturing the impacts of climate and land use changes on water supply [12][13][14].Analyses using global MSD models, including GCAM, typically rely on only a few alternative realizations of water supply (e.g.[6,[15][16][17][18][19]) due primarily to the lack of a large ensemble of GHM runs available.This lack of availability is possibly due to the computational burden of running and bias-correcting the global climate models that produce the inputs for hydrologic models [20][21][22][23].
A stochastic model of future runoff with global coverage can ameliorate a variety of existing challenges for the assessment of the multisector consequences of water availability.Stochastic watershed modeling is a recent approach that pairs deterministic watershed models with statistical models of their inputs or errors to generate many plausible synthetic streamflow realizations.These realizations can be used to help water managers account for uncertainty and risk in infrastructure planning and management [24].Stochastic watershed modeling has not yet been deployed globally.A global application of stochastic watershed modeling would allow for a probabilistic representation of future global water supply and for the identification of critical scenarios not represented by the small set of existing GHM realizations.In addition to the small number of GHM runs typically available for uncertainty propagation, a major challenge for hydrologic models-including GHMs-is that they tend to underestimate variance in streamflow and runoff, leading to misrepresentation of extremes (e.g.drought or flood events) [25].Stochastic watershed modeling addresses this issue by adding variability to the hydrologic model outputs and producing ensembles of streamflow or runoff with corrected statistical properties [24,25].
In this work, we apply stochastic watershed modeling to an existing GHM to generate alternative realizations of basin-scale runoff with global coverage.We use these realizations as inputs to GCAM, propagating uncertainty in water availability across the interconnected energy, water, and land use sectors.By linking many stochastic water futures to an MSD model, we explore the impact of hydrologic uncertainty and extremes on regional water and agricultural outcomes.

Methods
In this section, we describe the global stochastic watershed model that we develop, the data that we use in our model, and GCAM, the MSD model that we use to simulate the multisector impacts of future runoff variability.

Overview of stochastic watershed model
There are three major challenges associated with developing a global stochastic watershed model: (1) when applying a global stochastic model, we must account for spatial correlation across different water basins; (2) GHMs are computationally heavy, so our stochastic approach needs to be computationally efficient and not require many GHM simulations to generate ensemble members; and (3) there is no dataset of observed streamflow or runoff with a long enough historical record and sufficient spatial coverage to train the stochastic model.In this section, we describe how our stochastic watershed model addresses these three challenges.
Our stochastic watershed model is inspired by the approach of [26], which builds a time series model of the hydrologic model's errors, then uses bootstrap sampling from this error model to generate a stochastic daily streamflow ensemble.This procedure adds variability back to the hydrologic model output to reduce error and improve performance at extremes.By generating the stochastic ensemble using the output of a hydrologic model, we can account for anthropogenic influences on hydrology (i.e.climate and land use changes) not possible with direct simulation of a hydrologic time series [24].Here, we develop a statistical error model to generate global ensembles of annual basin runoff across 235 basins.We model runoff rather than streamflow because our MSD model (GCAM) uses annual runoff as the upper limit of renewable water supply.
Figure 1(a) illustrates the steps that our stochastic watershed model takes, which are described in detail in supplementary data S1.The model has three inputs: a historical reference, a model of annual runoff, and a projection of future annual runoff.The historical reference can be either observational data or a benchmark model that best approximates reality if insufficient observational data is available (see section 2.2).The model of annual runoff can be any deterministic hydrologic model with the appropriate spatial and temporal resolutions.Here, we use the Xanthos hydrologic emulator [27,28] as the model of future annual runoff (not as the benchmark) but another model could be used.We first calculate the errors between the historical reference and model for each basin and fit a linear regression to the errors with the model runoff as the independent variable.To generate the stochastic ensemble in the historical period, we sample errors from a multivariate normal distribution for each year and add the error back to the model.A multivariate normal distribution is used to account for spatial correlation of errors across the basins, thus addressing the first challenge mentioned above (see supplementary data S1 for explanation of construction of covariance matrix).To generate the stochastic ensemble for the future period, we use our linear error model to generate errors based on the future projection of annual runoff.Given its stochastic nature, our model does occasionally produce negative annual runoff values in low-runoff basins (0.32% of all flow values).We set these negative values to zero to reflect that there is negligible available runoff.GCAM operates in five-year time steps, so we apply a five-year backwards rolling mean to the annual runoff values which further reduces the impact of zero flows; across all flow values inputted into GCAM in this study there is only a single occurrence of zero flow (in the arid Sinai Peninsula basin in 2045).We assume that the spatial correlation of the errors across basins is the same as in the historical period.Our approach addresses the second challenge laid out above because it is computationally efficient and does not require additional runs or biascorrection of a GCM or GHM.

Runoff datasets
The third challenge that our stochastic watershed model confronts is that there is limited available observational data for annual runoff with global coverage.Therefore to train our error linear regression model, we use a well-known GHM, WaterGap2 [29][30][31], as the reference to which we compare the performance of the global hydrologic emulator Xanthos [27,28], the model that prepares water supply inputs for GCAM.Using a model as the reference when no observed dataset exists has precedent in global hydrologic and climate model comparisons [27,32].By using a model as our benchmark, we have a long record of data from 1901 to 2019.Future projection of annual runoff for 2020-2100 was produced using Xanthos forced with climate data from the CMIP6 GFDL global climate model with bias-corrected scenario SSP370 [21,33].All models output monthly runoff at a 0.5-degree grid resolution with global coverage, which we aggregate to annual runoff at the basin scale due to the coarser spatial and temporal resolutions of GCAM.Hereafter, we refer to WaterGap2 as the reference and Xanthos as the deterministic model.Additional detail on these GHMs is provided in supplementary data S2.These models are used in this study to introduce our new methods for exploring the multisector impact of hydrologic variability but alternative models or datasets could also be applied to our methods.
To evaluate the performances of the deterministic model and the stochastic watershed model in the historical period, we compare their performance to the reference using the Nash-Sutcliffe Efficiency (NSE) for each basin [34].NSE is a measure of model performance that captures the agreement of model predictions to observations, with an upper bound of one for a perfect model.An NSE of zero indicates that the model has no more predictive power than the reference's mean value.In addition to NSE, we consider the percent errors in mean and variance of the models relative to the reference.

MSD model
To explore the multisector impacts of uncertainty in future hydrology, we use GCAM, an integrated model of the socioeconomic, climate, energy, water, and land systems [8].GCAM is an open-source, dynamicrecursive partial equilibrium model that captures the long-term responses of the energy, water, and land use sectors to climate, socioeconomic, technological, and policy changes.The model operates in fiveyear time steps for 2020-2100 (with calibration years 1975-2015) and divides the world into 32 energyeconomy regions, 235 water basins, and 384 land-use regions (formed by intersecting the energy-economy regions with water basins).Each of these divisions is referred to as a market.GCAM solves for an equilibrium price in each market at each time step when supply equals demand for the energy, water, and land use sectors.The model uses a logit choice methodology to determine the shares of different energy and agricultural technologies in each market based on price and profit [8].
As demonstrated in figure 1(b), we use our stochastic runoff ensemble to modify GCAM's representation of water supply, which is divided into renewable water (surface water and renewable groundwater minus environmental flow requirements), nonrenewable groundwater, and desalinated seawater.Specifically, we take a five-year backwards rolling mean for each basin for each member of our stochastic runoff ensemble for 2020-2100 and use these values as each basin's maximum renewable surface water supply in GCAM.Cost curves for renewable water and non-renewable groundwater are used to define their costs; renewable water is the cheapest due to its relative ease of access while groundwater is more expensive due to extraction cost.Desalinated water has a fixed high cost and is used only for municipal and industrial purposes when other water supplies are limited.GCAM outputs a water price for each basin that represents the economic value of the last unit of water used to produce goods.While in reality water is often not traded in markets and does not have a price, the water price quantity is a useful measure for exploring important water sector dynamics such as preferential water allocation or increasing price that reflects water scarcity and allowing interaction with other sectors.GCAM has been used in previous studies to explore the physical and economic consequences of changing water supplies and demands under different socioeconomic, climate, hydrologic, and land use conditions (e.g.[4,6,35].Aside from the alternative water supply curves input to GCAM, this work follows the GCAM reference scenario which assumes moderate population growth and no major economic shifts or technological advancements compared to historical patterns [36]. Changing the renewable water supply input (i.e.propagating different realizations from the stochastic model) to GCAM will impact both the availability and the affordability of water resources for agricultural and energy uses.This can lead to changes in both the quantity of production as well as where production is occurring.For example, if there is less renewable water available in a basin in a given time period, that basin may need to tap into its available groundwater resources, which comes at a higher cost.As a result, the basin may shift to growing less irrigation-intensive crops.Since global demand for those crops remains relatively the same (we do not change any demand-related assumptions in this study), other basins in the world that have more renewable water available may shift to growing more irrigation-intensive crops.To quantify the multisector impacts of varying renewable water supply, we consider total water withdrawals, surface water withdrawals, water withdrawals for irrigation, water price, and irrigated cropland outputs from GCAM across multiple realizations.We quantify how these outputs are impacted by variability in water supply using a relative variability indicator that we calculate as the coefficient of variation in each output divided by the coefficient of variation in runoff.This relative variability ratio indicates a basin's ability to absorb the impacts of variability in runoff for each output.For example, a relative variability ratio greater than one for irrigated cropland for a selected basin would indicate that the uncertainty in irrigated cropland is amplified by runoff variability.A relative variability ratio less than one would indicate that the basin is able to absorb the impacts of runoff variability for irrigated cropland.

Deterministic model performance
We begin with an evaluation of the deterministic model (Xanthos), finding that it misrepresents variance compared to the reference in most basins.Figure 2(a) displays the percent error in mean and variance of the deterministic model compared to the reference for the 235 water basins considered in this study.In most basins, Xanthos has small error in the mean; 88% of basins have less than 20% error in the mean.In contrast, Xanthos struggles to capture the variance of the reference flows, returning greater than 20% error in half of the basins.These findings are consistent with those of [25] who considered much finer temporal and spatial resolution models.To illustrate the consequences of error in the model variance, we focus on the Indus and Mekong basins.Figures 2(b) and (c) show flow duration curves, commonly used in hydrology to illustrate the probability of exceeding runoff values of a given magnitude.These plots also report the NSE for the deterministic model relative to the reference.
Both the Indus and Mekong are examples of basins for which the deterministic model performs well on average but suffers at representing the extremes.Both basins have NSE values greater than zero (0.77 for the Indus and 0.46 for the Mekong, figures 2(b) and (c) and low error in mean behavior (1% and 2%, respectively, figure 2(a)).However, each basin's performance suffers at the extremes due to under-or over-estimation of variance.The Indus basin underestimates variance (−47%) while the Mekong basin overestimates variance (96%).When variance is underestimated, like for the Indus basin, the deterministic model tends to overestimate the lowest flows and underestimate the highest flows (figure 2(b)).Relying on the deterministic model in such cases would result in modeling the lowest flow years (droughts) and highest flow years as less severe than the reference indicates.Conversely, when variance is overestimated like for the Mekong basin (figure 2(c)), the model predicts that the lowest flow years are more severe than the reference and that there is more water available in high flow years than is actually the case.Both types of error in variance can lead to poor understanding of hydrologic extremes, which are linked to some of the most societally consequential outcomes across sectors [37][38][39].This motivates the need for stochastic models that can better represent uncertainty in the runoff prediction while also correcting for model bias and misrepresentation of variance.

Stochastic model performance in historical period
To evaluate the performance of the stochastic model, we begin by comparing an ensemble of 1000 stochastic realizations to the deterministic model and reference in the historical period.The ensemble contains runoff values for all 235 basins for all 119 years in the historical period.Figure 3(a) displays the median percent error in mean and variance across 1000 stochastic ensemble members, with both statistics showing substantial improvement across all basins compared to the deterministic model.The median of the stochastic ensemble has percent bias near zero for all basins, indicating that the stochastic ensemble is able to correct bias in the deterministic model (as shown in figure 2(a)).The stochastic ensemble also captures the variance of the reference better than the deterministic model, with 98% of basins having less than 5% error in median variance across the ensemble.In addition to average and extreme behavior, we also assess the models' ability to capture the behavior of the reference using Pearson's correlation coefficient.Figure 3(b) shows that the stochastic ensemble does not capture correlation as well as the deterministic model, which is expected because we are adding error to the deterministic model.We also find that the stochastic ensemble reasonably captures spatial correlation across the basins during the historical period; most differences in basin-to-basin correlation pairings fall within −0.2-0.2 and 70% of basin-to-basin correlation pairings have less than 10% error (see figure S3).See figure S4 for maps showing the performance of the deterministic and stochastic ensemble for all basins.
In figure 4, we highlight the Indus, and Mekong basins to show how the stochastic model corrects for misrepresentation of variance in the deterministic model.In terms of year-to-year variability, the stochastic ensemble captures the general pattern of variation in the historical reference while increasing the range of magnitudes at each time period by adding uncertainty (figures 4(a) and (b)).As indicated by figures 4(c) and (d), the stochastic ensemble also corrects for bias in the mean and variance in the deterministic model.The stochastic model is able to correct for bias in the mean and variance because we are including the error structure in the model prediction.It is important to highlight that the correction for bias in both mean and variance occurs whether the deterministic model is overestimating or underestimating these statistics compared to the reference (e.g.variance for Mekong and Indus basins).By correcting for bias in the variance, the stochastic ensemble better captures behavior at extremes than the deterministic model (figures 4(e) and (f)).

Stochastic model performance in future period
We can use this well-calibrated stochastic watershed model to explore variability in runoff under future  climate change.We apply our stochastic model to the future period 2020-2100 to create an ensemble of projected global hydrology at a basin level.Figure 5 shows a stochastic ensemble of future runoff for the Indus basin.Compared to the historical deterministic model (blue cross), the deterministic projection of future annual runoff (red cross) for this basin has lower mean and higher variance.The stochastic ensemble captures this shift, as reflected by the corresponding shifts in mean and variance of the future stochastic ensemble compared to the historical stochastic ensemble (figure 5(b)).In the historical period, the deterministic model for the Indus overestimates low runoff values and underestimates high runoff values; because we assume stationarity of errors, the model assumes this is the case for the future deterministic projection as well (figures 5(b) and (c)) and adjusts runoff values accordingly.

Multisector impacts of variability in future runoff
To assess the multisector impacts of uncertainty in future runoff, we run the stochastic ensemble 10 000 times and select 100 scenarios to run through GCAM.The scenarios are selected by uniformly sampling the distribution of cumulative runoff in the Indus basin from 2070-2100 to ensure that our selection covers the full range of the ensemble, including extremes (figure S6) since the goal of this exercise is an exploratory analysis.We focus on the Indus for sampling due to its unsustainable water use both historically and projected into the future [40][41][42].Each scenario includes the time series for the other 234 basins as well.In this section, we illustrate how our stochastic ensemble of water supply can be used to explore the impacts of alternative future runoff scenarios on water and agricultural resources.We begin   A basin with a relative variability ratio less than one is able to absorb the impacts of variability in runoff, a basin with a relative variability ratio of one has equal variability in withdrawals or irrigated cropland to variability in runoff, and a basin with a relative variability ratio greater than one experiences higher variability in withdrawals or irrigated cropland than in runoff.
with a global-scale assessment of basin-scale impacts globally followed by a discussion focused on the Indus.

Basin-scale impacts on global water and agricultural sectors
In figure 6, we show plots of relative variability, calculated as the coefficient of variation in total water withdrawals (figure 6(a)) and irrigated cropland (figure 6(b)) relative to the coefficient of variation in runoff at end-of-century (2100).The relative variability ratio indicates a basin's ability to absorb the of variability in runoff.For total water withdrawals, we find that all but one basin have a relative variability less than one, with many basins having values just above zero indicating that those basins have abundant water to meet demand, no matter how much runoff is available.Basins that do have a higher relative variability for total withdrawals (e.g.Sinai Peninsula = 1.09,Iran = 0.96, Northwest Mexico Coast = 0.90, Indus = 0.83) tend to be those that are already using all of their available water resources (both surface water and groundwater), so a difference in runoff more significantly impacts how much water they withdraw.Many of the basins with high relative variability are in countries associated with economic development driven by groundwaterintensive irrigated agriculture, such as Mexico and India [43][44][45].Many of the same basins with high relative variability ratios for total water withdrawals also have high relative variability for irrigated cropland (figure 6(b)), illustrating the importance of water availability for producing agriculture.There are 11 basins that have irrigated cropland relative variability ratios greater than one, indicating that uncertainty in runoff leads to greater uncertainty in land use.This finding emphasizes the utility of MSD models that represent the interactions between sectors when exploring hydrologic uncertainty.This analysis of relative variability in water demands due to variability in future runoff demonstrates how leveraging a large ensemble of water supply scenarios can be used to identify basins that are most susceptible to variability in future supply.We are able to highlight basins that may need to plan future infrastructure development to account for potential vulnerability of not being able to meet demands for water if their access to groundwater is unsustainable.

Multisector impacts of runoff variability in the Indus basin
We explore how our stochastic ensemble of future runoff can be used to characterize the basin-scale impacts of runoff variability on the water and agricultural sectors, using the Indus basin as an example.Figure 7(a) summarizes relative variability (defined as the coefficient of variation in a quantity of interest divided by coefficient of variation runoff, as explained in section 2.3) for 2020-2100 for four quantities: surface water withdrawals, water price, irrigated cropland allocation, and total water withdrawals for irrigation.We find that variability in surface water withdrawals is approximately one across time, indicating that surface water withdrawals are strongly tied to how much runoff is available.This reflects how the Indus basin struggles to meet its demand with available supply: when there is less runoff available, there are fewer surface water withdrawals because the basin is extracting as much surface water as it is able to at all times across all scenarios.Additional withdrawals come from more expensive sources such as groundwater, which is reflected in the high relative variability of the water price metric.In the Indus basin, water price varies across the years but is generally greater than one.Because water price reflects the price of the last unit of water used, its variable price indicates that the amount of more expensive groundwater this basin needs to use varies across the scenarios depending on how much runoff is available.For both irrigation water withdrawals and irrigated cropland, variability is always less than one but increases over time, especially after 2040 when groundwater in the basin is largely depleted across all scenarios (figure S7).Irrigation-related quantities become increasingly tied to available runoff as time goes on as groundwater resources in the basin are depleted [46].This finding highlights the benefits of combining stochastic watershed modeling with GCAM's detailed representation of both surface water and groundwater resources and their uses in other sectors.Without integrating these two modeling techniques, it would not be possible to assess the temporal evolution of the variability of the agricultural sector's response to runoff variability and the importance of groundwater in mediating the agricultural sector's response.
An additional benefit of integrating stochastic watershed modeling with MSD modeling is that we can identify scenarios of extreme conditions for a basin and explore the cascading multisector impacts of those conditions across space.As a demonstrative example for the Indus basin, we select the highest, median, and lowest runoff scenarios in 2100 for this basin across the ensemble of 100 scenarios that were run through GCAM (figure 7(b)) and explore how irrigated cropland allocation is impacted when the Indus experiences particularly high-or low-runoff years (figures 7(c) and (d)).Under the highest and lowest runoff scenarios for the Indus, the basin has a 7% increase and 8% decrease in irrigated cropland relative to the median scenario, respectively, reflecting the impacts of increased or decreased surface water on land use and crop production.In general, we find that the most impacted basins tend to be those that are also near the Indus geographically, although basins in North America and Africa have noticeable changes in percent of cropland allocation.While specific global land use dynamics cannot be definitively attributed to the runoff changes in the Indus basins, this finding nonetheless highlights the importance of multisectoral global modeling to explore how these cropland relationships emerge.By generating a large ensemble of scenarios, we show that we are able to explore the multisector impacts of extreme conditions that we cannot do if we only rely on a few marker projections.The analysis shown here for the Indus is just one example of how this work can be applied; future work can leverage our stochastic ensemble to select scenarios of interest to explore the possible multisector impacts of co-occurring extremes across multiple basins.

Conclusion
Available future water supply is expected to have a significant role in the evolution of the demands and supplies of land resources over the next century.There is a high degree of uncertainty in future water supply due to unpredictable climate and socioeconomic changes.Thus, modeling alternative water supply outcomes is critical for a better understanding of the multisector consequences of uncertainty in future water supply.While a large body of previous work has explored how different climate and socioeconomic assumptions impact water and agricultural availability, these studies tend to rely on only a few projections of future water availability.To our knowledge, there has never been a systematic study with global coverage that explores how variability in the hydrologic system's response to climate may impact the future supplies and demands of water and agricultural resources.
Here, we address this research gap by integrating stochastic watershed modeling and MSD modeling to generate a large ensemble of alternative scenarios of future runoff with global coverage and explore how variability across this ensemble is absorbed or amplified across the water and agricultural sectors.We introduce a novel stochastic watershed model for annual runoff that can operate across multiple basins and generate a user-specified number of stochastic realizations.Unlike more complicated approaches of stochastic modeling that rely on Monte Carlo sampling of the many parameters of a complex hydrologic systems model, our model uses a simple postprocessing approach that is based only on the errors between a historical reference and the deterministic model.We show that our model is able to capture mean, variance, and spatial correlation of the historical reference while correcting for any misrepresentation of average or extreme behavior that occurs in the deterministic model.The model is also easily applied to future projections of runoff assuming stationarity of model errors.
We demonstrate the utility of generating a large stochastic ensemble of future annual runoff by running an ensemble of 100 realizations through the GCAM MSD model.By combining these two modeling approaches, we can explore how variability in runoff impacts water and land use and identify the basins that are most or least robust to variation in runoff.We use the Indus basin as an example of a waterlimited basin where variations in surface water supply has increasing importance on cropland allocation and irrigation water demands as nonrenewable groundwater is depleted.We also show how GCAM's representation of global trade can yield additional insights about the multisector impacts of local droughts.
One limitation of the proposed stochastic watershed model is that it relies on the historical relationship between simulated runoff and error, and requires extrapolation beyond observed errors in some cases.A second limitation is the occasional generation of negative runoffs.These limitations could result in unrealistic runoff projections in some cases.While climate change is expected to increase the occurrence of shocks such as flooding or sub-annual droughts, we are unable to model the multisector impacts of these events due to the coarser temporal resolution of GCAM.We also note that input hydrologic data which global hydrology models rely on can be sparse and of low quality, so it is possible that the models we use are more accurate in data-heavy regions like the United States and Western Europe.Interpretation of the numerical results for decision making purposes should be done cautiously.In addition, MSD models like GCAM are limited in the interactions that they can represent and do not typically represent water rights.In the absence of that information, assumptions need to be made.One such assumption is preference for surface water before groundwater based on cost, which simplifies the complicated relationship between water rights and accessibility.
While our approach is applied here to annual runoff, it can be easily modified to generate a stochastic ensemble for streamflow or for coarser temporal resolutions.For finer temporal resolutions, an autoregressive component would likely need to be included due to seasonality or time lag in the historical errors.Future work could create more advanced approaches for dealing with model extrapolation and negative stochastic runoff values, such as the log-ratio approach developed by [26].There is also potential opportunity for combining the stochastic hydrology framework we employ here with explicit considerations of climate and socioeconomic uncertainties (i.e.model and scenario uncertainties as well as internal variability) to explore the interacting influences of these three drivers on future energy, water, and agricultural resources.

Figure 1 .
Figure 1.(a) is a schematic of stochastic watershed model developed for generating ensemble of future annual runoff across multiple basins.Red ovals indicate model inputs and blue boxes indicate calculations.(b) is a diagram that shows the integration of the stochastic watershed model with GCAM.Red boxes indicate model output (from stochastic watershed model or GCAM).

Figure 2 .
Figure 2. Performance of deterministic model (Xanthos) relative to reference (WaterGap2) in historical period 1901-2019.(a) is a scatterplot of mean (horizontal axis) and variance (vertical axis) in annual runoff across all basins for the deterministic model.(b) and (c) show flow duration curves for the Indus and Mekong basins, respectively.Nash-Sutcliffe efficiency (NSE) scores for the deterministic model are indicated on the flow duration curves.

Figure 3 .
Figure 3. Performance of 1000-member stochastic ensemble relative to deterministic model and reference in historical period 1901-2019.(a) is a scatterplot of mean (horizontal axis) and variance (vertical axis) in annual runoff across all basins for the stochastic ensemble.(b) shows kernel density estimations of linear correlation with the reference for the deterministic model and stochastic ensemble.Note: the axes of (a) are limited to highlight main results.Four out of 235 basins are excluded, but are shown in figure S5.

Figure 4 .
Figure 4. (a) and (b) Time series plots, (c) and (d) scatterplots of mean versus variance in annual runoff, and (e)-(f) flow duration curves for the Indus and Mekong basins for the historical period.The black and blue crosses in (c) and (d) represent the reference and deterministic model, respectively.Median Nash-Sutcliffe efficiency (NSE) scores for the 1000-member stochastic ensemble are indicated on the flow duration curves.

Figure 5 .
Figure 5. (a) Time series of annual runoff, (b) scatterplot of mean versus variance of annual runoff, and (c) flow duration curve for the Indus basin that demonstrate stochastic model performance for the future period for a 1000-member stochastic ensemble.(a) and (b) also show stochastic model performance for the historical period.

Figure 6 .
Figure 6.Maps of variability ratios for (a) total water withdrawals and (b) irrigated cropland relative to runoff in 2100.A basin with a relative variability ratio less than one is able to absorb the impacts of variability in runoff, a basin with a relative variability ratio of one has equal variability in withdrawals or irrigated cropland to variability in runoff, and a basin with a relative variability ratio greater than one experiences higher variability in withdrawals or irrigated cropland than in runoff.

Figure 7 .
Figure 7. (a) Time series of variability ratios for surface water withdrawals, shadow price of water, irrigated cropland allocation, and water withdrawals for irrigation.(b) Time series of stochastic runoff for Indus basin 2020-2100.Maps of percent change in cropland allocation compared to the median runoff scenario from (b) in all 235 basins in the (c) lowest and (d) highest runoff scenarios for the Indus basin.