Watershed monitoring always raises the same practical questions: “Where should sensors be placed?” “Which hillslope would be most adapted to perform a field experiment or deploy monitoring equipment?” or “Given a particular site, how representative is it of the watershed?” These are all challenging questions for critical zone and watershed sciences (Strachan et al., 2016). Watersheds consist of diverse biotic/abiotic above-/below-ground compartments interacting with each other through hydrological, ecological, geomorphological, and geochemical processes (Sivapalan, 2006; Wagener et al., 2007).

To answer this type of questions, there have been several studies on sensor placement optimization for soil moisture and snow distribution in watershed science (Kerkez et al., 2012; Oroza et al., 2016) and stream flow depths (Zhang et al., 2019). However, they have focused on estimating the spatial distributions of a single or two target variables. A recent pressing issue has been to observe and quantify interacting ecosystem or watershed responses to the disturbances resulting from climate change, such as coupled soil–plant responses to droughts or wildfires (e.g., Sloat et al., 2015; Mirus et al., 2017; Siirila-Woodburn et al., 2021; Newcomer et al., 2023), or to investigate those disturbances’ key modulators, such as subsurface conditions, which are often difficult to characterize over space and time.

There have been significant advances in quantifying the above/below-ground watershed compartments over space based on a suite of remote sensing datasets, including plant species distributions and traits (e.g., Madritch et al., 2014; Falco et al., 2019; Chadwick et al., 2020; Falco et al., 2020), soil thickness and soil properties (e.g., Patton et al., 2018; Yan et al., 2020), and bedrock variability (e.g., Parsekian et al., 2015; Uhlemann et al., 2022). The co-variability of these compartments has been documented based on analyzing multiple remote sensing and spatial data layers (e.g., Wainwright et al., 2015; Pelletier et al., 2018; Devadoss et al., 2020; Hermes et al., 2020; Enguehard et al., 2022; Wainwright et al., 2022a). In addition, dynamic properties (such as snow distribution and plant productivity) can be estimated across the watershed scale using airborne and satellite remote sensing (Painter et al., 2016; Dong et al., 2019).

In parallel, integrated ecohydrological and reactive transport models have been successfully implemented to describe and predict watershed behaviors from hillslope to continental scales (e.g., Bierkens et al., 2015; Carroll et al., 2019; Xu et al., 2022). Coupled three-dimensional surface water and groundwater flow models have been augmented by land models representing plant water uptake and evapotranspiration (e.g., Condon and Maxwell, 2015; Maina et al., 2021). These models can physically predict the water budget and export at the watershed scale. In addition, reactive transport models have been developed at the watershed scale to quantify geochemical processes such as weathering, nutrient cycling, and solute export (e.g., Li et al., 2017; Li, 2019; Maavara et al., 2021; Xu et al., 2022). These models can predict distributed properties such as soil moisture, water table depth, weathering rate, and other properties over space and time.

Despite the availability of spatial data layers, there have been few studies developing a systematic framework to optimize observation locations (including sensor/well locations) or to choose particular hillslopes/catchments to maximize the efficacy of critical zone experiments. In particular, there is a need to take advantage of landscape heterogeneity to identify key controls, capture distinct ecohydro-biogeochemical regimes, and/or test particular hypotheses with sites are properly distributed (e.g., Sloat et al., 2015; Sanders-DeMott et al., 2018). At the same time, a number of restrictions have to be considered, such as protected areas (e.g., private property and wilderness designated), road accessibility, disturbances to recreation, and others.

This study aims to develop an effective strategy for guiding the selection of observation/experiment locations by combining remote-sensing data layers and model simulation results. In particular, we explore how best to locate multiple observation sites across a watershed to investigate underlying mechanisms that dictate the sensitivity or vulnerability of ecosystems to climate disturbances such as droughts. The ecosystem and plant responses can be directly observed by satellite images such as the year-to-year variability of the Normalized Difference Vegetation Index (NDVI; e.g., Knowles et al., 2017, 2018; Dong et al., 2019; Wainwright et al., 2020). The underlying processes, such as soil moisture and groundwater table dynamics, cannot be observed over space, but we may be able to take advantage of model simulations.

We employ a two-step process: (1) regression analysis of satellite-derived ecosystem sensitivity as a function of model-derived hydrological state variables (soil moisture, water table elevation, snow depths, snowmelt timing, and melt rate) to identify key variables for the sensitivity and (2) cluster analysis based on these variables to identify regions that have divergent responses. This work leverages zonation studies by Wainwright et al. (2022a) and Maina et al. (2021). We consider the hillslope as a unit since the critical zone or watershed science studies are often based on the hillslope-scale observations/experiments, which capture the hydrological connectivity from the ridge to the valley.

Here, we study the East River and adjoining three watersheds near Crested Butte, Colorado, United States (Hubbard et al., 2018) as a demonstration site for the approach combining remote sensing and hydrological model data products to identify distinct hillslope zones that exhibit divergent eco-hydrological characteristics. With regard to the ecosystem response, we consider the heterogeneity of the foresummer drought sensitivity based on the Landsat NDVI developed by Wainwright et al. (2020). This drought sensitivity represents the plant productivity responses to early snowmelt and subsequent drought conditions in the early growing season. In addition, we use the hydrological state variables, including soil moisture and snow depth, computed by a hydrological model (Carroll et al., 2022b). The snow-soil–plant interactions have been studied by ground-based monitoring (e.g., Sloat et al., 2015; Sanders-DeMott et al., 2018) and a combination of ground-based and remote-sensing data (e.g., Devadoss et al., 2020; Hermes et al., 2020). We aim to improve our understanding of these interactions—particularly their spatiotemporal variability.

The proposed framework employs a two-step approach: (1) the regression analysis of the ecosystem sensitivity as a function of model-derived hydrological state variables, and (2) the clustering analysis to group the hillslopes with particular ecosystem functions. The spatial data layers are first spatially aggregated into hillslope-scale metrics. In Step 1, we also aim to improve our understanding of snow-soil–plant interactions over space and time.

In this study, we extended the approach by Wainwright et al. (2022a) and Maina et al. (2021) to include both key observed data from remote sensing on plant productivity (i.e., foresummer drought sensitivity) and model-computed hydrological variables (soil moisture and snow depth). We used the same clustering approach to compress multidimensional data and identify zones that capture co-varied properties. This approach enables us to explore the effect of water limitation on plant productivity as well as snow-soil–plant interactions over space and time. In addition, we added an extra step to identify the most representative and accessible hillslopes for future sampling to enable a model-experimental (ModEx) design [Department of Energy (DOE), 2023].

The results show that soil moisture variability is the best predictor of foresummer drought sensitivity. This is consistent with the current understanding of the water limitation on plant productivity in this region (Sloat et al., 2015). At the same time, this fact also adds more confidence to the model results. Although soil moisture and snow thickness are known to be highly heterogeneously impacted by microtopography (Falco et al., 2019; Devadoss et al., 2020), the model results capture the general trend and variability of soil moisture across the watershed. However, the floodplain area in the Slate watershed shows a surprising discrepancy, in the sense that while the soil moisture and elevation are both low, the sensitivity is not quite as high as in other regions with similar soil moisture values. Since its upstream hillslopes tend to have significant snow accumulation at high elevations, it would be possible that this low-elevation region benefited more from stream and/or subsurface water from the upper watersheds. This discrepancy suggests the need for further investigation into this region and potentially some model refinements.

The predictive performance of foresummer drought sensitivity is higher with hydrological variables (e.g., ΔVWC) compared to the static explanatory variables (such as elevation, geology, and other parameters used in Wainwright et al., 2020). This is potentially a result of the non-linearity among soil moisture, snowmelt timing, and elevation. The model-derived soil moisture physically aggregates multiple parameters and effects, such as elevation, aspects, and hydrological gradient, into a single variable that is more effective in describing the water limitation on plant productivity. Although the obtained zonation map is quite similar to the one based only on the remote sensing data in Wainwright et al. (2022a,b), there are some differences. This zonation map based on both remote sensing and modeling results has more zones at higher elevations (Zones 1–3). This is attributed to the non-linearity, such that the snow-free timing is more heterogeneous than other variables influenced by radiation at similar elevations. This also suggests that, while currently the elevation and elevation-based life zones (such as montane, subalpine, and alpine) are often used to represent the distinct ecosystem zones, model-based metrics can be used to further refine these distinct zones.

Additionally, the combination of satellite observation data and modeling results may provide alternative hypotheses for snow-soil–plant interactions. The snow-soil–plant interactions under climate change have been extensively studied in the recent decade (e.g., Schmidt and Lipson, 2004; Brooks et al., 2011; Sanders-DeMott et al., 2018). Different mechanisms and processes are considered, such as increased soil warming (Melillo et al., 2011), microbial or biogeochemical impacts (Schmidt and Lipson, 2004; Brooks et al., 2011), and the decrease in snowpack reducing insulation and increasing freezing-induced root damage (Sanders-DeMott et al., 2018). In our study domain, as mentioned above, plant productivity is considered to be limited by water availability, such that early snowmelt leads to early soil drying and creates drought conditions before the summer monsoon (Sloat et al., 2015; Dafflon et al., 2023). Although this study confirms this hypothesis, it also suggests a non-linear interaction across space and time and the presence of threshold values.

In our results, soil moisture is correlated with snow-free timing, but the relationship is non-linear. In the high-elevation north-facing region, as well as in the big snow year, the soil moisture does not depend on snow-free timing because soil moisture becomes saturated due to the deep snowpack, limited soil storage, and the lateral transfer of water from high alpine regions to upper subalpine zones (Carroll et al., 2019). In dry years and low-elevation regions, the correlation is more significant, suggesting these conditions and/or locations are soil moisture-limited. Based on these results, we may hypothesize that these regions, which are currently resilient to year-to-year climate variability, may become more sensitive to the variability in the future as the temperature increases and snow depth decreases at higher elevations. Since such non-linearity is difficult to detect through sparse measurements, spatially continuous modeling results are likely the most effective way to identify non-linearity and threshold processes.

We would note that, in this study, we did not consider the antecedent conditions from the previous year. We assumed that it could be limited in this system because, due to snow and thin soils in the montane-alpine region, the soil becomes fully saturated after snowmelt every year. Such a soil storage limitation has been reported in other regions (Hahm et al., 2019), and its impact on plant dynamics would be interesting to investigate in subsequent studies.

Our approach also identifies the different zones—i.e., the hillslope groups—that potentially have divergent snow-soil–plant interactions. In addition to the low- and high-sensitivity regions, clustering identified the intermediate- to higher-elevation regions where the topographic aspects are considered to play key controls. As global warming continues, snow-free timing is expected to decrease, but may or may not do so below the threshold in these regions since the topography is not likely to change in the same time frame as climate. In this future scenario, the steep hillslopes and the large variability of aspects and potential radiations would increase the spatial heterogeneity of drought sensitivities.

This type of monitoring design methodology is expected to be increasingly important as sensor and telecommunication technologies improve and a wide variety of fixed/permanent sensors are available (e.g., Hubbard et al., 2018). The U.S. Geological Survey’s Next Generation Water Observing System, for example, plans to have distributed sensors and observational sites for soil moisture, snow-water-equivalent, and others (US Geological Survey, 2023). In addition, sporadic ground-based surveys such as snow course measurements are still important to cover a larger space between fixed sensors and to create site-specific relationships with remote sensing data (Dressler et al., 2006). In particular, our approach is promising to guide the optimal deployment of wireless sensor networks (e.g., Oroza et al., 2016; Wielandt et al., 2023), in which the value of sensing various properties of interest might need to be balanced by the site accessibility, sensor resolution and accuracy, power requirement, wireless connectivity, total data cost, and complementarity to remote sensing products and other datasets.

Finally, this study shows the power of combining remote sensing data and modeling results to locate key observation locations and generate hypotheses that can be applied to other processes. Although plant productivity and trait data are increasingly available over time and space, subsurface properties are largely considered to be a “black box” (Sevanto et al., 2020). The modeling results are particularly powerful in this regard, providing subsurface state variables such as soil moisture. In addition, our method can determine the most representative hillslope in each zone that has road/trail access, which facilitates the selection of experiment/monitoring sites. The consideration of road/trail access could be easily extended to incorporate additional objectives in the site selection process. Now that ML approaches are widely available and more people are trained to use ML for extracting critical patterns in spatial data layers, model-guided monitoring design and hypothesis generation can be more readily employed in critical zone and watershed sciences.

Publicly available datasets were analyzed in this study. This data can be found at: the datasets are located on the U.S. Department of Energy repository ESS-DIVE: the simulation data link (https://data.ess-dive.lbl.gov/datasets/doi:10.15485/1889747) and the remote sensing data link (https://doi.org/10.15485/1841262).

HW performed statistical analyses and wrote the main text. NF processed remote sensing data layers and contributed to the manuscript. ES-W and RC provided model results and contributed to the manuscript. BD, IB, and YW provided perspectives on monitoring and sensor network and contributed to the manuscript. All authors contributed to the article and approved the submitted version.