Carbon emissions exacerbating the Western US water crisis

Water has long been a limited resource in the western US (WUS). In this semi-arid region, snowmelt and atmospheric evaporative demand act as the primary climatological controls on dry season water resources, the majority of which are used to irrigate a multi-billion-dollar agricultural industry^1,2,3. The bulk of precipitation for most of the region falls in the winter as snow^4,5. In the warmer months, snowmelt from mountainous regions contributes more than 50% of total surface runoff⁶, thereby acting as the primary source of surface water^7,8,9. Atmospheric evaporative demand and agricultural water use also peaks in these warm months, the latter of which accounts for about 80% of developed water use across the WUS^{2,10,11,12,13}. However, observed climate trends have strained water resources. In recent decades, reductions in mountain snowpack^7,14, earlier runoff timing^15,16,17, and increases in atmospheric evaporative demand^11,18 have increased the seasonal gap between local water supply and demand¹⁹, which is projected to widen further in the future⁹. While some of the changes in regional hydroclimate are due to natural variability, rising temperatures due to anthropogenic climate change (ACC) have already contributed to the region’s multi-year droughts^20,21,22, reduced spring snowpack²³, earlier snowmelt runoff^24,25, increased evaporative demand^21,22, and increased irrigation needs²⁶.

A complex water storage and conveyance system has been engineered to bridge the seasonal and geographic mismatch between water supply and water demand in the WUS¹. However, local demand still often exceeds available supply, leading to overallocation of surface flows and increased reliance on groundwater^27,28,29. Multiple factors contribute to groundwater overdraft. Reductions in surface water supply have accelerated groundwater pumping to satisfy irrigation demand, particularly during drought years, while simultaneously hindering groundwater recharge. Thus, recent decades have seen a rapid depletion of groundwater^30,31, raising the alarm about a coming water crisis^31,32. This is particularly evident in California’s Central Valley, which has experienced significant groundwater depletion over the past few decades³⁰. While this depletion is in part due to the long-term overallocation of groundwater, it has accelerated in recent decades, with depletion rates now 30–45% faster during the twenty-first century megadrought period compared to the late twentieth century^10,30. This has resulted in dry domestic groundwater wells and land subsidence³³, expensive legal challenges to surface water curtailments during dry years^27,34, and the passage of California’s Sustainable Groundwater Management Act (SGMA) with the goal of curbing groundwater use.

Although policies like SGMA target direct groundwater use, such efforts do not address the contributions of high emitting entities to climate change-driven water cycle disruptions. Thus, in response to emerging climate impacts, including this water crisis, several legal and policy mechanisms are being pursued to hold fossil fuel producers accountable for their contribution to such impacts and for their broad deception and disinformation campaigns^35,36,37,38, including ‘climate superfund’ legislation in California and county- and state-level lawsuits across the WUS^39,40,41. These legal and policy mechanisms rely on attribution methodologies to quantify the share of observed impacts attributable to emissions traceable to fossil fuel producers. Fossil fuel producers have contributed substantially to ACC: roughly 70% of global industrial carbon dioxide (CO₂) emissions since 1854 are traceable to 122 fossil fuel producers and cement manufacturers (“Carbon Majors”), and ~97% of these emissions occurred after 1950, a decade when many of these entities began to understand the consequences of their products^42,43,44. Recent studies have quantified the contribution of the Carbon Majors’ emissions to global temperature rise and sea level rise^44,45, ocean acidification⁴⁶, increased burned area in western North America’s forests⁴⁷, extreme heat events⁴⁸, and loss of gross domestic product⁴⁹. However, other direct impacts–including on regional water resources–remain unquantified. To fill this gap, there is a debate as to whether source-to-attribution studies are necessary (“attribution approach”), or if impacts can be estimated based on the proportion of total emissions that are traceable to the Carbon Majors (“proportional approach”)^49,50. This debate revolves around whether impacts can be scaled linearly using existing climate attribution analyses by a ratio of the increase in CO₂ emissions or Global Mean Temperature (ΔGMT) due to the Carbon Majors, or if these impacts are not linear enough to satisfy the assumptions of the proportional approach⁴⁹.

Here, we quantify impacts on WUS (west of 105°W) water resources–including water supply, water demand, and groundwater—attributable to the Carbon Majors, and compare the attribution and proportional approaches. We focus on three variables to assess WUS water: April 1st snow water equivalent (SWE_Apr1) and warm season streamflow (streamflow_warm) on the supply side, and annual total irrigation demand on the demand side. We focus on these months as April is near peak SWE for most basins, and hence is important for spring and summer streamflow, when local demand is at its highest (Fig. S1). We then assess how the climate change-forced signal in these variables propagates to groundwater decline.

To quantify the changes in these variables attributable to the Carbon Majors, we designed a source-to-impact attribution approach (see “Methods”). Leveraging estimates of ΔGMT forced by Carbon Majors emissions versus ACC from previous work⁴⁴, we first related ΔGMT to local meteorological variables using pattern scaling and global climate models (GCMs). Second, we used high-resolution meteorological data along with the pattern scaling results to develop counterfactuals which exclude changes forced from Carbon Majors emissions versus those from total ACC. Third, we translated the observed and counterfactual versions of these variables to SWE, streamflow, and irrigation demand using offline modeling. Finally, we used these variables to model changes in groundwater for a hotspot basin to further isolate the contributions of Carbon Majors and ACC on recent declines in groundwater. Aggregating across WUS basins, we estimate that Carbon Majors emissions since 1950 (CM1950) have driven declines of 15% in April 1 SWE and 6% in warm season streamflow, and a 2.4% increase in irrigation demand, accounting for 40–64% of the total ACC-driven changes across the region. These changes have widened the seasonal gap between water supply and demand, with decreased streamflow and increased irrigation demand during the warm season. For California’s agriculturally intensive Central Valley, our results suggest that the changes in water supply and demand due to CM1950 account for about a third of the observed, climate-driven groundwater loss during 2003–2024.

Contribution of carbon majors to changes in western water

Comparing hydrologic simulations of observed and counterfactual forcing (i.e., attribution approach), we show that emissions traceable to the Carbon Majors since 1950 (CM1950) led to a 15% (~35 km³ per year) reduction in April 1st SWE (SWE_Apr1) across the WUS over the past decade (2014–2024) (Fig. 1a, b), compared to a 36% (~83 km³ per year) reduction in SWE_Apr1 for all ACC (Fig. S2a). Thus, the reduction in SWE_Apr1 due to CM1950 accounts for 42% of the total reduction in SWE_Apr1 due to ACC, or 40–52% across each of the five macroscale (HUC2) basins. The largest total reductions in SWE_Apr1 are found in the Pacific Northwest (PNW) and California (CA) basins (Fig. 1a). Over the past decade (2014–2024), CM1950 has reduced SWE_Apr1 by 5.1 and 27.0 km³ in CA and PNW, respectively, each equivalent to 18–20% of the observed values for each basin (Fig. 1a, b).

We similarly found a 6% (~15 km³ per year) decline in warm season (April–September) streamflow (streamflow_warm) across the WUS over the past decade attributable to CM1950, compared to a 13% (~32 km³ per year) reduction attributable to ACC (Figs. 1c, d and S2b). The reduction in streamflow_warm due to CM1950 accounts for 47% of the reduction in streamflow_warm due to ACC, or 43–50% across each of the five HUC2 basins. As with SWE_Apr1, the largest absolute reductions in streamflow_warm are found for CA and PNW, with 4.0 km³ (13%) and 9.9 km³ (6%) less annual streamflow_warm, respectively (Fig. 1c). Notably, these decreases are smaller in magnitude than that of SWE_Apr1. While SWE is sensitive to temperature and precipitation⁵¹, streamflow reflects changes in total precipitation, evapotranspiration, energy balance, and snowmelt timing⁵². We also demonstrate a significant change in the timing of streamflow due to CM1950. The centroid of timing of runoff is on average ~5.6 days earlier across the WUS, and up to 30 days earlier in some mountainous regions such as the western slopes of the Washington Cascades (Fig. S3a). Correspondingly, we found an increase in cool season (October–March) streamflow (streamflow_cold) due to CM1950 for most basins (Fig. S3b), accompanying the forced decrease in streamflow_warm. The one exception is for the Upper Colorado Basin, which has experienced a decrease in both streamflow_cold and streamflow_warm. (Figs. 1b and S3b). Totaled across the two Colorado basins, the absolute decline in streamflow due to CM1950 is 3% (~0.43 km³ per year). While this is far lower than the declines seen for CA and the PNW, given the overallocation of the river, such small changes in streamflow can have outsized effects on the local water budget^28,53.

We found average annual increases in irrigation demand of 2.4% (~0.89 km³ per year) across the full study region over the past decade due to CM1950 (Fig. 2), compared to ACC-driven increases of 4.3% (~1.56 km³ per year; Fig. S2d). This irrigation demand is equivalent to the additional irrigation that would need to be applied to a crop in order to match the increase in evaporative demand over that field, assuming perfect irrigation efficiency. The increase attributable to CM1950 accounts for ~57% of the increase due to ACC (Fig. 2b), or 51–64% across each of the five HUC2 basins. While the absolute values are substantially smaller than the declines in SWE_Apr1 and streamflow_warm, they are no less concerning. Only a fraction of runoff reaches lower elevations due to evapotranspiration in mountainous regions, and of that, even less is stored and available for irrigation withdrawals during the warm season⁶. Furthermore, irrigation demand is at its peak in the warm months, coinciding with lowest runoff values (Fig. S1). Thus, an increase in irrigation demand in the warm months increases water stress during these normally dry months, further impacting agricultural water resource allocation.

Across basins, the largest absolute increases in irrigation demand due to CM1950 are found in the PNW (central Washington and southern Idaho) and CA (in the Central Valley of California) (Fig. 2a), corresponding to the most intensively irrigated agricultural regions in the WUS (Fig. 4a). Over the past decade, there has been ~0.26 km³ (2.8%) and ~0.36 km³ (1.9%) more annual irrigation demand due to CM1950 in PNW and CA, respectively (Fig. 2b). However, the strongest relative increase is seen over the PNW and Upper Colorado Basin (UCB) (Fig. 2b), coinciding with a larger increase in reference evapotranspiration (ETo) attributable to CM1950 (Fig. S4b). This finding suggests that areas with relatively lower existing irrigation demand may nonetheless face increasing biophysical water requirements due to rising evaporative demand, even without significant expansion of agriculture.

Proportional versus attribution approaches

We found evidence of nonlinearity in the changes in water supply and demand variables in response to ΔGMT, leading to differences between the proportional and attribution approaches (Fig. 3). For certain areas in the WUS, there is a notable difference between the two approaches for estimated changes in SWE_Apr1, streamflow_warm, and irrigation demand attributable to CM1950 (Fig. 3a–c), at times greater in magnitude than the uncertainty stemming from the ΔGMT and pattern scaling estimates in our attribution analysis (see “Methods”). For both SWE_Apr1 and streamflow_warm, these differences are a function of winter (December–March) temperature (Fig. 3g, h). For cold regions (winter Tavg < −8 °C), there is little to no difference, though also little attributable decrease (Fig. S5). However, for regions with near-freezing temperatures (0 °C <winter Tavg <4 °C), the proportional approach yields substantially larger decreases in SWE_Apr1 and streamflow_warm than the attribution approach (Fig. 3g, h). The rain-snow transition is more sensitive to warming, thus, the sensitivity of SWE to temperature is nonlinear^23,54. As the proportional approach cannot capture this nonlinearity, it overestimates the sensitivity of SWE_Apr1—and as a by-product, streamflow_warm—to temperature at these rain-snow transition elevations (Fig. S5). Conversely, the proportional approach consistently underestimates the increase in irrigation demand (Fig. 3c). As irrigation demand is a function of ETo, values increase nonlinearly with warming (due to the saturation vapor pressure curve). Therefore, each incremental increase in emissions (and thus change in temperature) results in larger increases in ETo. The proportional approach cannot capture this, instead assuming a linear relationship between temperature and ETo change, and thus underestimates irrigation demand, particularly for warmer areas (Figs. 3i and S5).

These differences remain when aggregating to HUC2 and HUC4 basins (Figs. S6–10). Totaled across the full study domain, the proportional approach underestimates irrigation demand by 0.14 km³ (19%) relative to the attribution approach, and overestimates SWE_Apr1 and streamflow_warm by 5.18 km³ (13%) and 0.58 km³ (4%), respectively (Fig. S7). For certain basins and time periods, these differences are greater in magnitude than the uncertainty estimates in our attribution analysis. For the CM1950 scenario, the difference between approaches in quantifying attributable changes in SWE_Apr1 and streamflow_warm for PNW and the full domain is within one measure of uncertainty stemming from the attribution analysis (10th–90th percentiles) but falls outside of the other measure (25th–75th percentile) (Fig. S7). Furthermore, these divergences compound when dealing with more recent time periods: for Carbon Majors emissions since 1990 (CM1990), all three variables (SWE_Apr1, streamflow_warm, and irrigation demand) fall near or outside of the 10th–90th percentile estimates from the attribution analysis (Fig. S8).

Combined effects

Hotspots for absolute changes in streamflow_warm attributable to CM1950 emerge in the PNW and the northern CA HUC4 basins (Fig. 4a), while irrigation demand hotspots are found in the Central Valley of CA, eastern Oregon and Washington, and southern Idaho (Fig. 4b). Four basins are ranked in the top 10th percentile for changes in both streamflow and irrigation demand: the Klamath, Sacramento, and San Joaquin basins in California and the Middle Columbia basin in Oregon and Washington (Fig. 4a, b). These basins have a large amount of irrigated agriculture and exist near the line separating water-limited versus energy-limited regions on the Budyko curve (Fig. S11).

For these basins, the changes in water supply and demand attributable to CM1950 have widened the window in which elevated demand cannot be met by surface water supply (Fig. 4c, d)^9,19. In these basins, irrigation demand has increased due to CM1950 during the warm months (May–October). Conversely, while streamflow has decreased during the warm season (June–September for the Middle Columbia basin, and April–July for the Sacramento basin), it has increased during the cool season. These seasonal changes have important implications for meeting water demand in the hot, dry months. Short of changing irrigation practices (e.g., deficit irrigation), shifting cropping practices (e.g., agricultural land transitions), or building out more storage, this generally forces an increased reliance on groundwater. Notably, the proportional approach is generally captured within the uncertainty bounds of estimates using the attribution approach, albeit with systematic underestimates of changes in irrigation demand and slight overestimates in summer streamflow changes.

To assess how these changes have affected groundwater in California’s Central Valley, we examined late-summer terrestrial water storage as a proxy for groundwater totaled across three HUC4 basins (Sacramento, San Joaquin, and Tulare basins). While we identified the Sacramento basin as in the top 10th percentile from both supply and demand terms (Fig. 4), because Tulare and San Joaquin are reliant on water imported from the Sacramento basin⁵⁵, these two basins are effectively also in the top 10th percentile for water supply. We fit a model to predict year-to-year changes in total water storage anomalies (and by proxy groundwater anomalies) as a function of streamflow_cold, streamflow_warm, and annual total irrigation demand (r² = 0.63; see “Methods”), reflecting that streamflow_cold and streamflow_warm contribute to groundwater recharge (and surface water availability), streamflow_warm can be used for irrigation (instead of groundwater), and irrigation demand can drive groundwater depletion (Figs. S12, S13).

In the Central Valley from 2003 to 2024, there has been an observed decline in groundwater of 33.9 km³. Assuming ~70% of this decline was due to overallocation, and ~30% due to climatic factors (inclusive of ACC and natural variability), observed groundwater depletion due to climatic factors was ~10.2 km³ during 2003-2024 (see “Methods”). We estimate cumulative declines of 6.5 km³ (10th–90th percentile estimates: 2.0–9.8 km³) and 3.5 km³ (1.9–4.7 km³) attributable to ACC and CM1950, respectively (Fig. 5a). Thus, ACC and CM1950 account for ~64% and ~34%, respectively, of the climatically-driven groundwater loss during 2003–2024, or roughly 20% and 10% of the total observed groundwater loss.

Notably, the proportional approach suggests lower cumulative decreases in groundwater (3.1 km³) due to CM1950 compared to the attribution approach (3.5 km³), but this discrepancy is small compared to the uncertainty in the groundwater reduction estimates (10th–90th percentile: 1.9–4.7 km³) (Fig. 5a, b). This is due to two factors. First, at the scale of the three Central Valley HUC4 basins, while the proportional approach underestimates the increase in irrigation demand, the change in streamflow_cold and streamflow_warm is similar between the two approaches (Fig. S9). Thus, the underestimation of groundwater declines in the proportional approach for this basin comes primarily from its underestimation of irrigation demand. Second, while there is relatively good agreement between the two approaches for CM1950, there is less agreement between the approaches for CM1990, suggesting the ability of the proportional approach to accurately estimate attributable changes varies with emissions time period (Figs. S10 and S14).

Emissions traced to the Carbon Majors have substantially reduced water supply and increased water demand across the WUS, with emissions since 1950 (CM1950) accounting for roughly half (~40–64%) of the impacts attributable to anthropogenic climate change (ACC) (Figs. 1 and 2). Specifically, CM1950 has led to ~35 km³/year (~28 MAF) less SWE_Apr1, ~15 km³/year (~12 MAF) less streamflow_warm, and ~0.9 km³/year (700,000 AF) more irrigation demand across the WUS over the past decade, corresponding to ~42%, 47%, and 56% of the changes attributable to ACC, respectively (Figs. 1 and 2). To put these values in context, the annual decline in SWE_Apr1 across the WUS due to CM1950 is roughly equivalent to the maximum capacity of Lake Mead, and the increase in irrigation demand due to CM1950 is ~1/3 of the annual residential water use of California. Our analysis revealed two spatiotemporal hotspots of hydrologic disruption concentrated in the Pacific Northwest and northern California, where basins have experienced a substantial decrease in SWE_Apr1 and streamflow_warm, and an increase in total annual irrigation demand (Fig. 4). In these basins, anthropogenic forcing has led to a widening seasonal gap between water supply and demand (Fig. 4c, d), which has profound impacts on groundwater. In the Central Valley, we estimated that ACC accounts for ~64% of the climatically-driven decline in groundwater from 2003 to 2024, while emissions from the CM1950 accounts for ~34% (Fig. 5), corresponding to at least ~20% and ~10% of total observed groundwater loss during this period. These results imply an increased reliance on groundwater, surface water storage, and conveyance systems to capture surface-flows when they are available, thus foreshadowing persistent issues with regional water availability and storage as temperatures continue to rise. While these changes largely mirror prior attribution studies that have considered total anthropogenic climate change^23,24,26, we add specific value in isolating changes attributable to Carbon Majors.

While our procedure quantifies uncertainty, it may not capture the full range of structural uncertainty arising from the selection of climate models, choices in parameterization, and hydrologic model structure. First, modeled temperature and precipitation co-vary, in that models with a warming bias may tend toward wetter or drier biases in different regions⁵². Consequently, our approach of quantifying climate model uncertainty through symmetrical percentile matching (i.e., cold bias matching with dry bias) will not capture compounding extreme alignments (i.e., hot-dry and cold-wet) that may be prevalent in parts of the WUS. Therefore, we are confident in our attribution analysis but remain agnostic to these potential biases in our uncertainty calculations, which is an area for improvement in future studies. Second, while there is uncertainty stemming from the choice of hydrologic model (VIC), this structural uncertainty should be second order compared to the climate change effects⁵⁴. Third, our groundwater counterfactual uses a simple, linear model, which may contain inaccuracies in how it resolves sensitivities between input variables (e.g., streamflow, irrigation demand) and groundwater anomalies. Furthermore, as this model does not account for water management choices (e.g., managed aquifer recharge, changing irrigation practices⁵⁶), it does not distinguish how much of the decline is due to environmental factors (such as climate) versus historical overallocation. In accounting for this, by apportioning 30% of recent depletion as due to climatic factors, we provide a conservative estimate of climate change-forced declines. An area of future research is to use a more sophisticated offline model to assess climate change-forced changes in groundwater. Therefore, while we believe that our procedure still provides a reasonable estimate of uncertainty, these uncertainty bounds can be improved with future research.

The effects of changing water supply and increasing water demand extend beyond agriculture, initiating a cascade of impacts across the WUS with implications for wildfire, ecosystem health, and the overall resilience of regional water systems and associated sectors. Increased atmospheric evaporative demand enhanced aridity and contributed to the ongoing megadrought across the southern WUS^21,57. Similarly, changing supply and demand terms also lead to higher wildfire risk in forests of the WUS, with increased evaporative demand driving increases in burned area^47,58, and decreased SWE and earlier snowmelt linked to larger fires in montane regions⁵⁹. Furthermore, decreased water supply in the warm season and associated increased reliance on groundwater can lead to groundwater wells going dry, particularly for residential wells³³. Moreover, earlier snowmelt and peak runoff can lead to larger magnitudes of peak runoff and contribute to flooding risk for downstream communities^60,61. Furthermore, even small absolute declines in the water balance can have serious implications for local water budgets. For example, in the overallocated Colorado River Basin, low precipitation and warm temperatures have resulted in historically low flows, leading to reallocation of water to states^28,53,62. While we find relatively smaller decreases in streamflow (largely comparable to previous modeling efforts^53,62) and increases in irrigation demand for this basin compared to other basins in the WUS, even a small shift in water availability and demand can have outsized effects on the regional water budget^28,53.

There are several schools of thought regarding whether contributions to localized impacts can be determined simply as a proportion of ΔGMT, or whether nonlinearities inherent to certain types of impacts require direct impact attribution^49,50. As the proportional approach leverages results from the attribution approach (Eq. 13), our analysis effectively tested how well the proportional approach does in replicating the results from the attribution approach. We found systematic over- and under-estimations in localized supply and demand variables from the proportional approach (Fig. 3). At elevations near the rain-snow transition, the proportional approach overestimates reductions in supply terms (SWE_Apr1 and streamflow_warm) compared to the attribution approach. Conversely, the proportional approach underestimates irrigation demand, especially for areas with high baseline values. Yet when estimating basin-wide changes (Central Valley groundwater), we found relatively good convergence for CM1950, with the proportional approach falling within the middle 50% of estimates from the attribution approach for CM1950 (Fig. 5a), though falling outside of those uncertainty estimates for CM1990 (Fig. S14a), suggesting that the proportional approach becomes less accurate as one looks at more recent time periods of emissions. Whether these differences are substantial enough to warrant impact attribution in lieu of proportional approaches depends on the application of these findings. Different disciplines (e.g., science, management, policy, law) require different levels of precision⁶³. The proportional approach estimates may be suitable for demonstrating causality but may be less suitable for damage estimation. An area for future research is to assess what types of impacts might respond more linearly to ΔGMT (thus lending themselves to the proportional approach) versus which are substantially nonlinear enough to require a source-to-impact attribution analysis, and how these uncertainties translate into policy and legal spheres.

Our results have direct implications for climate governance and accountability. By quantifying the contributions of major emitters to regional hydrologic stress, this study provides both methodological advances and scientific evidence to emerging liability cases and climate ‘superfund’ legislation^40,41. As many individual entities within the Carbon Majors had knowledge of the projected climate effects of their products as early as the 1950s⁴³, the detrimental impacts to WUS water availability that we identify here highlight the consequences of the Carbon Majors’ continued fossil fuel and cement production in recent decades. By quantifying the role of major industrial emitters in worsening regional water stress, this study offers further evidence to inform economic damage attribution, climate risk disclosure, and future liability frameworks.

To estimate the change in water resources (water supply, water demand, and groundwater) attributable to all anthropogenic climate change (ACC) versus emissions traceable to the Carbon Majors, a source-to-impact attribution approach was designed (Fig. 6). First, we leveraged prior work using a reduced complexity climate model to determine the increase in global mean temperature (ΔGMT) attributable to four forcing scenarios (ACC, and Carbon Majors emissions after 1854, 1950, and 1990) (Step 1). Next, pattern scaling was employed to calculate gridded, monthly relationships between four meteorological variables (reference evapotranspiration, or ETo; minimum and maximum temperature, or Tmin and Tmax; precipitation, or Pr) and ΔGMT (Step 2). Counterfactual estimates of ETo, Tmin, Tmax, and Pr were then calculated using the ΔGMT estimates from Step 1 and pattern scaling results from Step 2 (Step 3). We used counterfactual climate data to model counterfactual snow water equivalent (SWE) and streamflow (Step 4) and crop irrigation demand (Step 5). Finally, groundwater was modeled as a function of these variables and actual and counterfactual groundwater estimated (Step 6).

Data

Our study region was set as the Continental US west of 105°W. For the study region, observational monthly minimum and maximum temperature (Tmax, Tmax), precipitation (Pr), relative humidity (RH), downward shortwave radiation (Rs), and wind speed (u₂) were retrieved from gridMET⁶⁴ for 1980–2024 to compute reference evapotranspiration (ETo). Additionally, daily Tmin, Tmax, u₂, and Pr were retrieved from the extended Livneh dataset for 1975–2018 for input into the hydrologic model^65,66. As this dataset is only available through 2018, we extended the data through 2024 using climatically aided interpolation⁶⁷ that leveraged nClimGrid-Daily⁶⁸ anomalies post-2018. This approach superimposes anomalies of daily Tmin, Tmax, and Pr based on a period common to both datasets (1981–2010) from nClimGrid-Daily to Livneh. Our extension of this dataset was shown to reproduce the observed variability over time (Fig. S15).

Monthly T, Pr, Td, Rs, and u₂ were retrieved from 114 simulations across 12 GCMs participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) from the historical (1850–2014) and SSP2-4.5 (2015–2100) experiments⁶⁹ (Table S2). The SSP2-4.5 scenario was used as it represents a “middle-of-the-road” scenario for emissions; moreover, as this is an attribution exercise focused on retrospective changes and scenarios do not significantly diverge until mid-century, suggesting the choice of future climate forcing scenarios has little bearing on our results. While there are different approaches to model selection⁷⁰, we chose these models because they have all the variables in the two experiments, have at least 3 ensemble members available, together have a mean equilibrium climate sensitivity of ~4.3 (Table S2), within the likely range⁷¹. These data were resampled to a common 1-degree grid. GMT was calculated for each of the 12 models from 1850–2100 by taking an area-weighted mean (at their native grid) of 2 m global surface air temperature.

Monthly total, CO₂-adjusted ETo was calculated using both gridMET and CMIP6 data following the reference crop Penman-Monteith equation, adjusted so that stomatal conductance decreases with increased CO₂ (Eq. 1)^72,73.

Step 1: Increase in GMT attributable to Carbon Majors

We draw from modeling⁴⁴, which used the MAGICC7 reduced complexity climate model to assess the impacts of emissions traced to the Carbon Majors on the increase in GMT (ΔGMT) for a series of historical counterfactual scenarios (Fig. 6, Step 1). MAGICC7 is a reduced complexity climate-carbon cycle model with a hemispheric upwelling diffusion ocean model which has been calibrated against observations and CMIP6 GCMs and was used to quantify temperature responses to emissions scenarios in the IPCC AR6 cycle^74,75,76,77. The simulations conducted in Sadai et al.⁴⁴ used the 2024 updated “Carbon Majors” database which includes scope 1 and scope 3 emissions traced to the largest 122 oil, gas, coal, and cement producing entities, including publicly-traded corporations and state-owned operations⁴².

We retrieved a probabilistic ensemble of 600 model realizations of ΔGMT for 1750–2024 for each of four scenarios: full anthropogenic forcing (ACC), and three counterfactual scenarios with emissions from the Carbon Majors removed from full historical emissions forcing (or from SSP2-4.5 emissions after 2015) during different time periods: 1854–2020 (CM1854), 1950-2020 (CM1950), and 1990–2020 (CM1990)^44,78. The scenario designated CM1854 reflects the full period of record, CM1950 captures the years after which it has been shown that the industry had knowledge of the climate effects of their products⁴³ and CM1990 reflects the years since the start of coordinated international efforts to address climate change⁴⁴. Here, we use the retrieved simulations^44,78 to recalculate ΔGMT for the mean as well as the 10th, 25th, 75th, and 90th percentile values (Table S1 and Fig. S16).

Step 2: Pattern scaling

The relationships between ΔGMT and the four variables (ETo, Pr, Tmin, and Tmax) were derived via pattern scaling using GCM-derived data at the common 1-degree grid (Fig. 6, Step 2). First, an average over all ensembles for each GCM was taken to reduce the influence of internal variability and better reflect the signal of forced change^79,80. Next, these ensemble averages were temporally smoothed using 31-year rolling means for each month (ETo, Pr, Tmin, and Tmax) or year (ΔGMT). This process further diminishes any signal of internal decadal variability that may remain after averaging across ensemble members to better elucidate the anthropogenic signal²³.

For each model (mod) and month (m), two calculations were run pixel-wise between the variable (ETo, Pr, Tmin, Tmax) and ΔGMT: Pearson’s correlation (1901–2050) to test for linearity (Fig. S17), and a linear regression to determine the change in each variable (v) per 1 °C increase in GMT (Eq. 2 and Fig. S18). The regression was fit on 1901–2050 for Tmin, Tmax, and Pr. However, for ETo, it was limited to 1980–2050 due to the aerosol effect that dominates the signal for ETo due to global dimming before 1980⁸¹. For Tmin and Tmax, the regression coefficient was used as the pattern scaling value, while for ETo and Pr, percent change was derived from the regression coefficient to use as the pattern scaling value. Finally, an additional, annual pattern scaling value was derived for each model using annual total ETo and Pr and annual mean Tmin and Tmax using the same approach (Eq. 3).

Previous work has indicated that pattern scaling results in acceptable estimates of anthropogenically-driven changes in VPD for the Western US (WUS)⁴⁷, ETo for California’s Central Valley⁸², extreme heat⁴⁹, and precipitation globally⁸³. However, to test whether the assumption of linearity between GMT and each variable holds across our study region, a second verification step was taken. Based on an assessment of the correlation maps (Fig. S17), time series of each variable were extracted for four counties that fell in the lowest correlation grid cells (Fig. S19). For these four counties, model-based values for each variable were extracted, spatial means (Tmin, Tmax) or totals (ETo, Pr) taken, and annual means (Tmin, Tmax, ΔGMT) or totals (ETo, Pr) derived and then smoothed using 31-year centered rolling means. From these, for each of the 12 model ensembles and for the multi-model mean, the raw change (Tmin, Tmax, GMT) and percent change (ETo, Pr) were calculated for each year (Fig. S19). Notably, our results supported the assumption of linearity, with one caveat. While the gridded pattern scaling estimates of ETo, Tmin, and Tmax all suggested a linear relationship, there were a number of grid cells for Pr with a lower correlation coefficient (<0.9; Fig. S17). However, these grid cells are in places with limited precipitation change (Figs. S17–S18).

Step 3: Counterfactual variables and uncertainty quantification

Counterfactuals of ETo, Pr, Tmin, and Tmax were derived by removing anthropogenically-forced signals from the observational data (Fig. 6, Step 3; Eqs. 4–8). Thus, the ACC counterfactual refers to all anthropogenic forcing removed (inclusive of Carbon Majors), while the CM1854, CM1950, and CM1990 counterfactuals refer to all emissions associated with the Carbon Majors removed after 1854, 1950, and 1990 (Table S1; Fig. S16). The counterfactuals account for uncertainty as a function of the probabilistic ensemble from MAGICC7 (600 realizations) and the 12 pattern scaling values (Figs. S20–S21).

For each variable (v, e.g., Tmax) and scenario (scen, e.g., CM1950), monthly multi-model mean pattern scaling values (B1) were multiplied by yearly means of each realization (r) of ΔGMT to derive a single, mean estimate of the change in the variable (Δv) due to the scenario for each month and year (Eq. 4). Next, annual pattern scaling values (B2) from the 12 models (mod) were multiplied by the 600 realizations (r) of ΔGMT for each year (y), yielding 7200 estimates of Δv for annual values (Eq. 5). From these 7200 estimates, the mean and the 10th, 25th, 75th, and 90th percentiles (pctl) were retrieved. Next, the difference between the mean and each of the percentile estimates for Δv.annual were additively removed from the multi-model mean Δv to derive 10th, 25th, 75th, and 90th percentile estimates for each month (Eq. 6). These resulting Δv values were bilinearly interpolated to the same resolution as the observational data (gridMET or Livneh). Finally, counterfactual values (v.cf) were calculated by removing Δv from the observed values, done additively for Tmin and Tmax (Eq. 7) and multiplicatively for ETo and Pr (Eq. 8). For the gridMET values, the monthly delta was removed from each corresponding month of observational data. For the Livneh values, the monthly delta was removed from each daily value of the corresponding month in the observational data.

Step 4: Water supply

To quantify forced impacts on the supply side of water availability, the Variable Infiltration Capacity model (VIC, version 4.2) was run for the WUS to SWE and streamflow (Fig. 6, Step 4). VIC is a well-accepted hydrologic model, and has been used in a number of hydrologic studies in the WUS^6,7,52,54,84. VIC was parameterized using the Livneh dataset at 1/16° resolution and all simulations were performed in a full energy mode at 3-h timestep with water years 1975–1979 as a model warm-up and 1980–2024 as a simulation period.

VIC was first run using daily observational Tmin, Tmax, u₂, and Pr in order to evaluate accuracy in reproducing observed trends over the study domain. The model was calibrated using the Livneh dataset and validated using observations of soil moisture, snow water equivalent (SWE), actual evapotranspiration (AET), and streamflow. Readers can refer to Livneh et al. ⁶⁵ for more information on model performance. We note that the choice of the hydrologic model can influence results. However, we anticipate this structural uncertainty to be second order compared to the climate change effects since our analyses are focused on snowpack and seasonal volumes, as opposed to low flows⁵⁴.

VIC tends to create artificial glaciers in some pixels because if temperatures are cool enough, a given year’s accumulated snow does not fully melt, leading to a perennial accumulating snowpack in the model. We opted to mask out these pixels so that artificial glaciers would not skew the results given their exceedingly high SWE values⁸⁵. These pixels have a negligible contribution to total streamflow volumes in the region, with even basin-level, summer contributions being small⁸⁶. Climatological maximum SWE is <4000 mm across SNOTEL stations and basins across the WUS⁸⁵. Therefore, we created a spatially fixed glacier mask to screen for potential artificial glaciers by setting pixels that ever reach >4000 mm to NA.

In addition to the observational run, VIC was run 20 times using combinations of counterfactual Tmin, Tmax, and Pr: for each of the 4 counterfactual scenarios (ACC, CM1854, CM1950, CM1990; Table S1) and each of the 5 uncertainty parameters (mean and 10th, 25th, 75th, and 95th percentiles). Monthly SWE and total streamflow (baseflow plus runoff) were retrieved from these runs, and from the monthly values, annual values derived: April 1 SWE; annual centroid of timing of runoff (CT), or the day of the water year (i.e., October 1-September 30) when runoff peaks²⁴; cold season total streamflow (October–March), and warm season total streamflow (April–September). These variables were then aggregated as volumetric totals to HUC2 and HUC4 basins for later analysis.

Step 5: Water demand

Next, irrigation demand across the WUS was calculated using ETo and Pr from gridMET (Fig. 6, Step 5). Assuming full irrigation to a field to meet transpiration demands, irrigation demand is estimated as the difference between crop water demand and available moisture, where crop water demand is derived as the monthly (m) product of ETo and an area-weighted crop coefficient (Kcw)^11,26,87 (Eq. 10). While our calculations assumed perfect irrigation efficiency, drip irrigation has an 85% efficiency with flood irrigation at an even lower efficiency.

A correction factor (D) was applied to counterfactual estimates of ETo (Eqs. 10 and 11). Aerosol emissions peaked in 1980, before a detectable influence of GMT on ETo: therefore, in 1980, the anthropogenic influence on ETo consisted of reduced shortwave radiation (from aerosols), increased CO₂, and no detectable change in temperature, resulting in a net negative influence on ETo until the turn of the century. While other approaches (e.g., a low-pass filter approach) can account for this effect, the pattern scaling reliance on GMT that we employ here is unable to capture the changes in anthropogenic forcing before 1980. We correct for this by taking the difference (D) between ETo_cf derived using pattern scaling (PS) versus a low-pass filter(LP) in 2000 (Eq. 11; Fig. S22). D was then scaled by the difference between each Carbon Majors counterfactual scenario to derive deltas to apply to those.

Crop coefficient (Kc) values were retrieved to create the Kcw layer⁷². While efforts have been made to update Kc values to reflect improvements in cropping practices (related to drought-tolerant cultivars, fertilization, etc), many of these efforts are done for relatively small regions and have not been robustly estimated across different crop types^87,88. Hence, the FAO Kc values, used in conjunction with CO₂-adjusted ETo, remain a credible and reproducible approach for use in larger, regional studies. From the FAO Kc values, monthly, area-weighted Kc values (Kcw) were calculated for each gridMET pixel for the WUS (Text S1; Table S3). This procedure uses Kc values and development length stages⁷², irrigation data (30 m) for 2017⁸⁹ to mask out non-irrigated areas, growing season calendars from USDA^90,91,92, and crop-specific land cover data (30 m) from the Cropland Data Layer (CDL) for a reference year of 2017⁹³.

The “available moisture” term in Eq. 10 was calculated using a simple monthly time step bucket model⁹⁴, with soil water holding capacity of 150 mm set for perennial crops and 50 mm set annual crops and forage⁸². This model accounts for fluxes of moisture into (precipitation) and out of (evapotranspiration) the soil with monthly available moisture being a combination of precipitation infiltration and plant available soil moisture. A 1-year warm-up using ETo and precipitation climatologies (1980–2024) was used to set up the storage term values.

Irrigation demand was calculated pixel-wise at the resolution of gridMET. The gridded CO₂-adjusted ETo and counterfactual ETo were multiplied by the monthly, gridded Kcw values, and the resulting crop evapotranspiration (ETc) values fed into the bucket model along with precipitation, to estimate irrigation demand under actual conditions versus counterfactual conditions (Eq. 10). In total, 21 estimates were created: observed values, and each of the 4 counterfactual scenarios (ACC, CM1854, CM1950, CM1990; Table S1) and each of the 5 uncertainty parameters (mean, and 10th, 25th, 75th, and 95th percentiles).

Step 6: Groundwater impacts

Groundwater was modeled as a function of water supply and demand terms. Monthly terrestrial water storage (TWS) was retrieved from GRACE-FO for 2002–2024⁹⁵. TWS was cropped to the three HUC4 basins that comprise the Central Valley of California (Sacramento, San Joaquin, and Tulare), and volumetric totals were taken over the basins. While TWS includes soil moisture, snowpack, reservoir and other surface water storage in addition to groundwater, much of the long-term change in TWS has been attributed to declines in groundwater in the Central Valley. Further in these basins, the vast majority of TWS at the end of summer is groundwater as snowpack is nominal and surface reservoir and soil moisture is at its annual minimum. Thus, we used August-September mean TWS values as a proxy for groundwater anomalies (Fig. S23). However, there are gaps in the GRACE product, including a 1-year gap corresponding to the transition from GRACE to GRACE-FO. To fill these gaps, a reconstructed GRACE product (GRACE-REC) was retrieved⁹⁶. Values in GRACE-REC were bias-corrected to GRACE + GRACE-FO by aligning the means (Fig. S23). Next, for each missing year, the percent change in the bias-corrected GRACE-REC data from the previous year to the missing year was derived and used to estimate the missing August-September value based on the prior year in the GRACE + GRACE-FO data.

While chronic groundwater overdraft due to overallocation has occurred for multiple decades, groundwater declines have accelerated since 2000, with these changes largely associated with climate-driven factors including significant multi-year droughts. Fitting trends on the groundwater declines from 1960 to 1999 versus 2000 to 2021 from Liu et al. ³⁰, we produce a conservative estimate of ~30% of groundwater loss since 2000 as climatically driven. Thus, groundwater anomalies were reduced to 30% of total as an estimate of climatically-driven groundwater declines. From these, the year-on-year change in groundwater anomalies was derived. A model was fit to the year-on-year change in groundwater anomalies using an ordinary least squares regression as a function of cold season (October–March) streamflow, warm season (April-September) streamflow, and annual total irrigation demand (adjusted r² = 0.63). Counterfactual groundwater anomalies were estimated using the coefficients and counterfactual irrigation demand and streamflow (Eq. 12).

To assess uncertainty in each counterfactual scenario (ACC, CM1854, CM1950, and CM1990), the four percentile estimates from ΔGMT and pattern scaling (n = 4) used to calculate counterfactual groundwater. Cumulative changes in groundwater were derived from the year-on-year changes, then counterfactual estimates were subtracted from the actual estimate, resulting in cumulative changes in groundwater attributable to each scenario.

Proportional versus attribution approaches

We compared two approaches to calculating the change in each variable due to Carbon Majors emissions: the “attribution” approach (comprising steps 1–6 in Fig. 6, and methods described above) versus the “proportional” approach. The proportional approach assumes that the changes in each variable scale linearly with ΔGMT, and thus the change in each variable attributable to the Carbon Majors can be derived based on the ratio of ΔGMT due to ACC versus due to Carbon Majors emissions (Eq. 13).

To derive the ratio between ΔGMT attributable to the Carbon Majors and ΔGMT attributable to ACC, the MAGICC7 results were subset to 2014–2024 and the mean ΔGMT retrieved for each scenario (Fig. S16). To estimate the change in each variable due to the Carbon Majors under the proportional approach, the change in the variable attributable to ACC was then multiplied by the ΔGMT ratio.