Resolution matters when modeling climate change in headwaters of the Colorado River

Water supplies in the Colorado River Basin (CRB) serve ∼40 million people, irrigate over 5.5 million acres, and drive turbines for more than 4200 megawatts of electrical generating capacity (US Department of the Interior 2012). Headwater catchments, though accounting for only 15% of the 640 000 km² area, provide ∼85% of CRB streamflow (Christensen and Lettenmaier 2007) (figure 1(a)). The volume of streamflow produced in these headwater catchments, coupled with faster warming at high-elevations (Pepin et al 2015), suggests that local-scale processes are of critical importance to determining future streamflow. These catchments are topographically complex, leading to steep temperature gradients and nonlinear relationships between water and energy fluxes. Regional and global scale models reduce topographic gradients (figure S1 (available online at stacks.iop.org/ERL/15/104031/mmedia)). Here, we use a localized, physically-based model that integrates groundwater, surface water, and land surface processes (Maxwell and Miller 2005, Kuffour et al 2020) to compare a suite of projected climate impacts as grid resolution is coarsened from 100 m to 1 km in a representative headwater catchment. We find that high-resolution models predict that 4 °C of warming increases headwater evapotranspiration (ET) 10% more than coarse-resolution models, thus predicting streamflow reductions that are 4% greater than coarse resolution models of the same catchment. This discrepancy is driven by increased sensitivity in high-elevation evergreen forests to warming temperatures, a local-scale process that is corroborated by observations (Goulden and Bales 2014) and missed at coarse modeling scales. These results suggest that streamflow predicted by coarse-resolution models may underestimate future, climate-induced reductions to downstream water supplies.

Zoom In Zoom Out Reset image size

Most water management decisions assume statistical characteristics of the historical record can be used to predict the likelihood of future conditions. This assumption of stationarity is not always valid, especially as climate change alters the underlying statistics of the system (Milly et al 2008). Assumptions of stationarity are often chosen as the best available option to address a highly uncertain system; however, when model physics are robust and explicit, a nonstationary model can be effectively used to understand nonstationary processes such as climate change (Serinaldi and Kilsby 2015). As a result, models that use the physics of the system instead of historical data trends have been applied at multiple scales in hopes of improved estimates of conditions in a non-stationary future climate. A remaining challenge is that there is significant uncertainty between models. For example, a review of future Colorado River streamflow studies predicts a range from −6% to −45% of current flows by mid-century (Vano et al 2014) across different modeling methodologies and platforms.

Recent discussions of model formulations address the role of increasing both process complexity and model resolution simultaneously. Global Climate Models (GCMs) have the advantage of incorporating broader climate patterns and land-atmosphere interactions. A particular challenge for GCMs is resolving complex topography. For example, at a resolution of 200 km the topographic gradients of mountain regions are less resolved than at higher resolutions (figure S1). This introduces difficulties for GCM's in high elevation regions (Beniston 2003, Leung and Qian 2003, Shen et al 2018, Polo et al 2020). While the spatial resolution of GCMs is improving to 25 and 50 km resolutions in some cases (Pascale et al 2017, Vecchi et al 2019), even at these resolutions resolving complex topography remains a challenge (figure S1). Moving from the global to local scale, models increasingly improve the representation of complex topography and processes driven by local heterogeneity in vegetation, soil, and geology. Water management models resolve more local processes compared with GCMs and RCMs.While these models improve hydrologic process representation, most are applied at regional-scale resolutions (1–50 km) and simplify physical processes with conceptual storage elements and routed overland flow so that large suites of scenarios can be run in short periods of time for operational use (Christensen and Lettenmaier 2007, Hrachowitz and Clark 2017). Despite calls for better topographic representations and higher resolutions (Wood et al 2011, Bierkens et al Bierkens, 2014), no high resolution, integrated modeling simulations of future CRB streamflow have been conducted in Rocky Mountain headwaters. Here, we extend the regional representation of headwater catchments to the local scale with 1 km and 100 m lateral resolutions to better understand the role of complex topography in altering streamflow predictions under climate change. Specifically, this study addresses whether modeling resolution affects predicted streamflow exports to the CRB from headwater catchments under climate change.

A headwater catchment that is representative of typical CRB headwaters is selected for this experiment given that it is not computationally feasible to run these simulations over all headwater catchments across multiple resolutions and climate scenarios. The model domain covers the East River basin in Colorado, USA. The 255 km² catchment exports water to the Gunnison River, which in turn contributes 40% of streamflow to the Colorado River at the CO-UT border (Spahr et al 2000). The East River is a high-elevation, mountainous, snowmelt-dominated headwater that is considered broadly representative of other headwaters in the Upper CRB given its geomorphic and land cover characteristics, bedrock composition, and diversity of stream orders and energetics (Battaglin et al 2011 Battaglin et al 2012, Pribulick et al 2016, Winnick et al 2017, Carroll et al 2018, Hubbard et al 2018). With elevations ranging from 2705 m up to 4123 m leading to 1418 m of relief, the study area encompasses multiple ecosystems- alpine, subalpine, montane, and riparian- and multiple hydrologic settings from high energy mountain streams to low-energy meandering floodplains.

The geology is dominated by Cretaceous age sediments of the Mancos Shale and Mesa Verde formations with Paleozoic and Mesozoic age strata—Morrison, Maroon, and Gothic formations—dominating its eastern boundary; igneous intrusions are present as laccoliths and stocks throughout the basin (Gaskill et al 1991). Landcover is dominated by evergreen forests including various species of spruce, fir, and pine, and high-elevation grasslands with intermingled aspen forests, shrublands, and rocky areas devoid of vegetation (Homer et al 2015) (figure 2). With the exception of Precambrian crystalline rocks and significant stands of lodgepole pine (Pinus contorta), the vast number of bedrock and land cover types found throughout headwaters of the Upper CRB are present within the East River basin again indicating its broad representativeness.

Zoom In Zoom Out Reset image size

From October to May most precipitation falls as snow in the East River, and from July through early September precipitation usually falls in heavy monsoon bursts. The domain encompasses three meteorological stations. At the high elevation, 3261 m, snow telemetry site (Schofield, site number 737), average annual precipitation is 1200 ± 230 mm (30-year period of record) and average annual temperature is 0.6 °C. At the mid elevation, 3097 m, snow telemetry site (Butte, site number 380), average annual precipitation is 670 ± 120 mm (30-year period of record) and average annual temperature is 2.1 °C. At the relatively low elevation, 2915 m, EPA Castnet site (Gothic, site number GTH161), annual precipitation is 640 ± 100 mm and average temperature is 1.8 °C (WY2006), though this site uses a tipping bucket gauge to measure precipitation, which generally underestimate snowfall. Both SNOTEL sites (Butte and Schofield) use separate transducers for snow water equivalent and precipitation, making these estimates more reliable than those at the Castnet station.

Two models of the domain were constructed at 1 km and 100 m lateral resolutions (Foster and Maxwell 2018). Both models have five subsurface layers discretized into 0.1 m, 0.3 m, 0.6 m (soil), 8 m (geology), and 21 m (weathered bedrock). The choice of five layers was based on balancing higher resolution near the surface for soil moisture relationships with vegetation, and computational expense in the 100 m model. This discretization is consistent with the Noah land surface model implementation (Chen and Dudhia 2001). The catchment covers 1420 m of vertical gradient from 2700 m up to 4120 m; at 100 m lateral resolution this gradient is reduced to 1390 m and at 1 km to 1200 m. At 100 m resolution, vegetation patterns with elevation are well captured, but at 1 km resolution these relationships are less clear (figures 2(a) and (c)). Similarly, the geology that is well represented at 100 m resolution, with regions of preferential groundwater flow defined by thin layers of sandstone and limestone, is simplified after aggregation to coarser resolution (figures 2(b) and (d)). Model development, including an extensive sensitivity analysis of hydraulic conductivity (K) and Manning's n for the same baseline climate year simulated here (WY2006), is documented in Foster et al (2018).

These two resolutions not only represent a move from a high-resolution model at a regional scale, 1 km, to high-resolution at the catchment scale, 100 m (figure 2), but aggregation to the 1 km model has a dramatic effect on inputs that control local hydrology—a near 200 m loss in elevation gradient, vegetation distribution with elevation is simplified, a loss of preferential groundwater flow paths, and a simplification of the stream network.

The physically-based ParFlow-CLM (Maxwell and Miller 2005, Kollet and Maxwell 2006) model was used in this study. Like other integrated hydrologic models (Camporese et al 2010, Therrien et al 2012), ParFlow solves for overland flow simultaneously to subsurface flow in both the saturated and unsaturated zones. ParFlow is fully coupled to CLM in order to completely represent snow processes, vegetation dynamics, and a full energy balance at the surface. The version of CLM implemented in ParFlow has been updated independently as new information on land surface schemes has become available (github.com/parflow). For example, changes to the ET parameterizations with respect to global ET partitioning (Lawrence et al 2007) have been implemented in ParFlow-CLM. Furthermore, transpiration, snow physics and soil evaporation partitioning have undergone sensitivity analysis and been validated against observations and to other models of regional to global scale partitioning (Jefferson and Maxwell 2015, Maxwell and Condon 2016, Jefferson et al 2017, Ryken et al 2020).

The model is driven by observed meteorological datasets from water year 2006 (SI 2, figure S2). Water year 2006 was chosen to represent baseline conditions because it most closely matched historical means of temperature and precipitation data available from the two SNOTEL sites within the domain (1985–2014) and the main purpose of the experiment is to explore the impact of resolution. Future work will explore dry years, wet years, and other conditions outside of the average. Twenty-six plausible climate change scenarios are developed by perturbing the atmospheric conditions from this baseline (SI 2, figure S2) according to likely scenarios for the Rocky Mountains in the IPCC report (Lukas et al 2014). Precipitation predictions are seasonally variable, so we test eight likely precipitation shifts ranging from −5% of summer precipitation to +5% summer and +10% winter precipitation. Each of these precipitation scenarios is run at three temperature conditions- 0 °C change (current conditions), +2 °C (average RCP4.5 projection for Colorado Rocky Mountains by midcentury), and +4 °C (maximum RCP4.5 projection by midcentury) (table 1). While other atmospheric variables will also be altered in a future climate, using perturbations isolates the impact of each (temperature vs. precipitation). Individual perturbations allow for a detailed tracing of the mechanism driving differences in modeled predictions. Studies have shown that elevated CO₂ in a future climate may increase water use efficiency, thus partly compensating for higher ET at higher temperatures (Guerrieri et al 2019, Reitz and Sanford 2019). Future work could explore the impact of model resolution using downscaled future climate predictions or explicitly addressing the role of increased CO₂, which would also provide opportunities to constrain the uncertainty in estimates of reduced streamflow presented here.

Each simulation was run repeatedly for multiple years using the same WY2006 forcing dataset (an average of 3 years/simulation) until reaching dynamic equilibrium, a state where the mass and energy balance did not change year-to-year. The 100 m model was run with 64 processors and required approximately 14 h per year simulated, the 1 km model was run on one processor and required approximately 4 h per year simulated.

Temperature changes have a larger impact on total water content than precipitation changes for all results (SI 3, figure S3). While temperature changes remain the dominant driver of water storage reductions, precipitation impacts increase in subsurface variables (figure S3: 4,5). While this is true in the average water year case, the impact of temperature vs. precipitation could change in wetter or drier years as the water availability shifts in the catchment, though Rocky Mountains west of the continental divide have been shown in previous research to be more energy- than water- limited (Foster et al 2016).

Baseline case streamflow shows a different pattern at different resolution for peak runoff (figure S3: 3), likely due to different snowpack development and melt patterns as wider elevation ranges are resolved in the 100 m model. The majority of elevation loss between the 1 km and 100 m resolutions occurs at high elevations that also remain much colder, pushing the melt timing to later in the year. Two peak runoff periods are observed in both the 1 km and 100 m model, but the 100 m model has a much larger second peak due to this effect of delayed snowmelt. While streamflow timing is critical to predicting flood risks and maximizing reservoir storage, the CRB is unique in having reservoir capacities that are over four times the average annual inflow rates (Vano et al 2014). Many management decisions can therefore be made on an annual basis. Resolution comparisons of climate scenarios on an annual scale demonstrate that the 100 m model is much more sensitive to temperature perturbations and similarly sensitive to precipitation perturbations than the 1 km model (figure 3).

Zoom In Zoom Out Reset image size

ET increases at a rate of 12% per degree of temperature increase in the 100 m model, while only increasing at a rate of 9% per degree in the 1 km model (figure 3(a)). These rates yield a 10% difference between the high and coarse resolution models possible by 2050. In average rates of ET across the domain, the increase is 29 mm yr⁻¹ per degree of temperature increase in the 100 m model, and only 22 mm yr⁻¹ per degree of temperature in the 1 km model (figure S4(c)). In the high-resolution model, precipitation impacts to ET increase at higher temperatures, indicating that as temperatures increase these headwater systems become more arid and more sensitive to shifts in precipitation. Differences in ET alter the water balance, resulting in changes to streamflow. At 4 °C of warming, the 100 m model predicts streamflow declines of 16% (from 2.6 to 2.2 cms, figure S4(a)), while the 1 km model predicts declines of only 12% (from 2.8 to 2.5 cms).

A major difference between the 100 m and 1 km modeling resolutions is the increase in potential for lateral groundwater flow at higher resolution. To test whether lateral flow accounted for reduced streamflow at higher resolution, we ran additional simulations after modifying the model physics to only resolve the vertical (not the lateral) flow of groundwater (SI 5). If, in fact, lateral flow were the mechanism, the divergence between resolutions should disappear in the altered physics case. Instead, we find that ET diverges by resolution at a faster rate when only the vertical component of subsurface flow is included. At coarse resolution, ET increases 11% per degree of warming and at high resolution ET increases 16% per degree (figure S5). This ET divergence with resolution as temperatures increase indicates that lateral groundwater flow mitigates water loss to the atmosphere in a warming climate (there is less divergence in the full physics case) and is, therefore, not the mechanism for reduced streamflow.

The 100 m model covers an elevation range that is 190 m larger than the 1 km model, and 184 m of this additional relief is at higher elevation than in the 1 km model. While these high elevations do produce slightly colder temperatures (figure S2(d)), the ET change over 4 °C of warming is similar between modeling resolutions (figure 4) at the upper limit of elevations (above 4000 m and generally above treeline). The major source of difference between the 100 m and 1 km models is the high elevation evergreen trees that are resolved at 100 m (figures 4(a) and 2(c)) and missed at coarse resolution (figures 4(b) and 2(a)). While the percentage of total evergreen forests change <1% due to resolution (figure 2 table), there are fewer high elevation evergreens due to topographic aggregation in the 1 km model. Transpiration and canopy evaporation measurements of evergreen forests in headwater catchments have been shown to be very sensitive to changes in temperature, with their ET signals driving reductions in headwater streamflow volumes in the Sierra Nevada (Goulden and Bales 2014). Similar results have been shown in the Alps as well, where a high-resolution ecohydrological model (250 m) demonstrated that increases in high elevation forest ET drives runoff reductions missed at coarse scales (Mastrotheodoros et al 2020). Here, ET from the high elevation evergreen trees is limited by temperature in the Baseline case, but at +2 °C and +4 °C ET increases nonlinearly compared to ET increases in lower elevation trees (figure 4(a)). After 4 °C of warming, the change in ET from the baseline model is much larger for evergreen forests between ∼3300 m and 3600 m than for other landcover types, including low elevation evergreen trees (figure 4). This impact is driven mainly by an increase in canopy sublimation during fall and winter months (SI 6, figure S7). It should be noted that there are large uncertainties in modeling snow sublimation, especially given the dearth of observations at these high altitudes. However, the results found here are consistent with high-elevation, sublimation-specific research (Gascoin et al 2013, Sexstone et al 2018, Ryken et al 2020). The 1 km model does not resolve the complex landcover patterns with elevation that drive this behavior (figure 2), resulting in a reduced average change in ET at coarse resolution (figure 3).

Zoom In Zoom Out Reset image size

Despite improvements in computational power and hydrologic models, it is still challenging to do a CRB-wide study at the 100 m resolution presented here, especially to incorporate extensive climate simulations. However, further exploration of the mechanism identified here could ultimately lead to an upscaling relationship, where this effect could be accounted for, even in lower resolution models. Previous work has demonstrated an upscaling relationship for hydraulic conductivity in topographically complex regions, and it is possible that a similar equation could be derived for landcover and ET relationships with topography (Foster and Maxwell 2018).

In summary, this study extends existing model predictions of climate change impacts on the CRB with a local-scale, integrated hydrologic model of a headwater catchment (figures 1 and 2). We show that evergreen forest sensitivities to temperature changes are nonlinear with elevation (figure 4), affecting estimated ET rates and streamflow export to downstream basins (figure 3). While the entire CRB has been shown to be sensitive to both precipitation and temperature (Hoerling et al 2009, Vano et al 2014), we find that headwater hydrology in the Rocky Mountains is more sensitive to increases in temperature than to precipitation changes predicted for the region in the next 50 years (figure S3). Given that 85% of Colorado River streamflow is generated in headwater regions, mostly on the west slope of the Colorado Rocky Mountains, water managers need localized and physically-based model predictions to accurately plan for the coming century as Southwestern populations increase. Furthermore, if these results are generalizable to other headwater basins, they demonstrate that regional and larger-scale models of the CRB may underestimate climate change impacts to headwater hydrology, thus overestimating future water supplies for the Colorado River.

This study was funded by a National Science Foundation Blue Waters Fellowship under award numbers ACI-0725070 and ACI-1238993, Department of Energy Office of Science Graduate Student Research Award under award DE-SC0014664, and as part of the Watershed Function Scientific Focus Area funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Number DE-AC02-05CH11231.

The data that support the findings of this study are available upon reasonable request from the authors.

LMF, KHW, and RMM conceived this study; LMF acquired funding for the first two grants listed, and KHW and RMM contributed to acquiring funding for the last grant listed; LMF conducted modeling, post processing, visualization, and data analysis while RMM supervised and advised; LMF and RMM interpreted results; LMF wrote the original draft and KHW and RMM reviewed and edited subsequent drafts.