Empirical attribution of a drying Himalayan river through remote sensing and secondary data

Streamflow regimes are rapidly changing in many regions of the world. The ability to attribute these changes to specific hydrological processes and their underlying climatic and anthropogenic drivers is essential to formulate effective water policy. Traditional approaches to hydrologic attribution rely on the ability to infer hydrological processes through the development of catchment-scale hydrological models. However, such approaches are challenging to implement in practice. In particular, models have difficulty capturing hydrological regime shifts, where changes in the dominant hydrological pro5 cesses alters the relationship among hydrological fluxes. Additionally, observational uncertainties might preclude closure of the catchment-scale water balance, which is a pre-requisite for most catchment-scale hydrological models. Here we present an alternative approach to hydrological attribution that leverages the method of multiple hypotheses. We generate and empirically evaluate a series of alternative and complementary hypotheses that pertain to hydrological change. These hypotheses concern distinct components of the water balance and are evaluated independently. This process allows a holistic understanding of 10 watershed-scale processes to be developed, even if the catchment-scale water balance remains open. We apply the approach to understand changes in the Upper Jhelum river, an important tributary headwaters of the Indus basin, where streamflow has declined dramatically since 2000 and has yet to be adequately attributed to its corresponding drivers. Using remote sensing and secondary data collected from the watershed, we explore changes in climate, surface water, and groundwater. The evidence reveals that climate, rather than land use, had a considerably stronger influence on reductions in streamflow, both through 15 reduced precipitation and increased evapotranspiration.


Catchment layout
The upper Jhelum river serves as one of the six main tributaries of the Indus river (Fig. 1a), supporting residents in both India and Pakistan. The watershed covers approximately 13,000 km 2 and is bounded by the Himalayan range to the northeast and Pir Panjal range to the southwest, both of which drain into the Kashmir valley. The river flows east to west along the Kashmir valley, passing through Wular lake before discharging at the outlet (Fig. 1b). The valley rests upon a layer of alluvial 135 sediment which partially overlays a layer of glacial till, the combination of which extends nearly 1 km into the subsurface (Dar et al., 2014). Both mountain ranges on either side consist of limestone Karst formations, and streamflow is supported both by groundwater recharge in the valley and karst springs emerging from mountainsides (Jeelani, 2008). The main stem of the upper Jhelum is intercepted by Wular lake (Fig. 1d) between gauges (iii) and (iv). A number of other, smaller lakes intersect tributaries within the watershed and seasonally inundated wetlands occupy much of the center of the valley (Fig. 1e). Agricultural land comprises the majority of the valley and supports rice paddy, maize, and wheat, as well as 145 noticeably increasing fruit orchards (DES, 2015). In addition to the above metrics for annual streamflow at the watershed outlet, streamflow observations were obtained from the The Department of Irrigation and Flood Control for Jammu and Kashmir at four additional stations (Fig. 1b) for the period 1955-2013. Standard data integrity checks (Searcy and Hardison, 1960) and flow separation (Nathan and McMahon, 1990) procedures were applied to ensure temporal and spatial consistency (see Supplementary Information). We used these data in 155 our evaluation of hypotheses, as described below. iv.
i. In these images (d-g), Landsat bands swir2-swir1-red are mapped to red-green-blue, so that water and snow are clearly visible as dark and blue pixels, respectively. MODIS imagery in (a) (Vermote, 2015) and Landsat composite imagery in (b, d-g) from U.S. Geological Survey were prepared and downloaded using © Google Earth Engine (Gorelick et al., 2017).  Figure 2. Hypothesis generation and evaluation. We develop a general empirical approach to attribution that (a) utilizes the water balance to (b) generate overarching research questions and (c) construct multiple hypotheses. The hypotheses are addressed (d) using remote sensing and in situ secondary data and (e) by applying empirical analyses, the results of which are (f) used to infer relationships regarding changing hydrological processes. We implement this approach in the Upper Jhelum watershed by defining the long-term water balance in terms of streamflow (Q), precipitation (P ), evapotranspiration (ET ), and changes in catchment water storage (∆S). Specific research questions (red), data (blue), and analyses (green) focus on the detection and characterization of changing water balance fluxes in the Upper Jhelum. The inference step is contingent on understanding how each flux has changed and provides a deeper understanding of watershed hydrological processes.

Implementing the method of multiple hypotheses
We use the water balance as a guiding framework to attribute the decrease in Jhelum streamflow. In particular, it serves as a conceptual tool to build alternative hypotheses regarding changes in each of the water balance fluxes that could explain the 160 observed reduction in streamflow (Fig. 2). Testing these hypotheses individually allows us to build understanding about the different pathways of hydrologic change throughout the catchment. Importantly, this approach does not rely on the predictive ability of an aggregate hydrological model, where calibration would depend upon on specific information that is not necessarily available. The multiple hypotheses are developed for each component of the water balance in the following paragraphs, along with approaches and datasets (also see Fig. 2, Table 1) to evaluate these hypotheses.

Hypotheses
Changes in precipitation (P ) could occur through different mechanisms that would have distinct impacts on streamflow. For instance, a reduction in precipitation would decrease the water balance inputs and reduce the water available to generate streamflow. Alternatively, an increase in storm frequency and reduction in average storm size could produce an increase in 170 vadose zone water storage, greater evapotranspiration (ET), and reduced streamflow (Zhao et al., 2019). We therefore pose the following hypotheses: -Hypothesis 1: Changing climate led to a reduction of annual precipitation.
-Hypothesis 2: Changing climate led to a change in rainfall seasonality.
The first hypothesis represents a direct reduction in the water input to the catchment that would generate a reduction in stream-175 flow. The second hypotheses represents a shift in the availability of water throughout the year. Such a shift in precipitation from a season with low ET to high ET could decrease average streamflow.
-Hypothesis 3: Changing climate led to greater storm frequency.
A reduction in storm size would reduce the "fast" component of streamflow (i.e., quickflow, McCaig, 1983) and combined with an increase in storm frequency the fraction of rain that is captured within the vadose zone would increase, leading to 180 greater bare soil evaporation and plant water use (Zhao et al., 2019). These changes would reduce quickflow and groundwater recharge, and thus ultimately reduce annual streamflow volume.

Data sources
To address these hypotheses, monthly precipitation data were obtained from the Indian Meteorological Department (IMD, Srinagar) for the period 1984-2013 at Srinagar, Kupwara, and Pahalgam stations (stations ii, vi, and vii in Fig. 1b). We 185 interpolated monthly precipitation to the watershed scale using Thiessen polygons, which we found were able to close the water balance better than less parsimonious approaches that account for elevation gradients (see Sect. S2, Fig. S3).
The monthly frequency of the IMD precipitation dataset precluded analysis of storm recurrence. Therefore, we utilized daily precipitation records from the gridded PERSIANN Climate Data Record (CDR, Ashouri et al., 2015). The purpose of PERSIANN CDR is to provide consistent precipitation data for long-term climate analysis dating back to 1984 at 0.25 degree 190 resolution. We used PERSIANN data as an indication of temporal variation in the distribution of storm size and frequency. In particular, we counted the number of storms each year in each pixel to see if storm recurrence had changed.

Hypotheses
Evapotranspiration (ET) affects streamflow by reducing the volume of water stored in the vadose zone that could have otherwise 195 produced streamflow. A reduction in streamflow could therefore be generated by an increase in ET, either through changes in potential evapotranspiration or vegetation properties and land use. We therefore include the following hypotheses: -Hypothesis 4: Climate change and warmer air temperatures led to greater evapotranspiration.
-Hypothesis 5: Land use change toward water-intensive crops led to greater evapotranspiration.
Both hypotheses are grounded in observed, ongoing changes within the watershed. Temperatures have been warming (Zaz 200 et al., 2019) and there has been a notable shift towards orchard plantations in portions of the valley (Romshoo and Rashid, 2014), which may use more water than traditional crops due to a longer growing season (Allen et al., 1998) and better access to subsurface water storage (Zhang et al., 2018). Both changes might have led to increased evapotranspiration and a reduction in streamflow.

205
We required an approach to estimate seasonal and annual evapotranspiration for the periods before and after 2000. Multiple approaches have been developed to estimate evapotranspiration (ET) using remote sensing products, but few provide robust estimates that would allow consistent comparisons of the the pre-and post-2000 periods. For instance, MODIS provides an 8-day ET product (Mu et al., 2013), but this dataset is only available since the launch of the Terra satellite in late 1999. The surface energy balance approach (SEBAL, Bastiaanssen et al., 1998) can be used to estimate ET from Landsat imagery prior 210 to 2000, but relatively few Landsat images were available during this time period and the instantaneous nature of SEBAL may therefore present biased estimates of seasonal ET.
We therefore applied a regression framework based on the concept of reference evapotranspiration (Allen et al., 1998), where evapotranspiration is defined as ET = ET 0 × k. In this equation, actual evapotranspiration (ET ) is characterized by a reference evapotranspiration (ET 0 ) and mediated by a crop coefficient (k), which respectively capture the effect of climate 215 and land use (i.e., vegetation). For calculating reference ET, we used the Hargreaves equation (Hargreaves and Samani, 1985), which provides estimates of ET 0 at monthly timescales or larger. The Hargreaves equation captures the effect of climate drivers and accounts for extraterrestrial radiation (R a ), air temperature (T A ), and the diurnal temperature range (T R ): Extraterrestrial radiation R a varies seasonally with changing solar declination angles which are associated with latitude and topography. The diurnal temperature range T R depends on a variety of climate conditions including humidity, soil moisture, 220 precipitation, and cloud cover (Dai et al., 1999;Geerts, 2003). We group these parameters into a single seasonal parameter, a s ∝ R a T 0.5 R , which represents mean climate conditions within each season and captures the seasonal variability in the Upper Jhelum watershed. We assume that the greenness of the vegetation canopy exerts the primary control on the stomatal conductance, such that we can define the crop coefficient as k ≈ N DV I c , where c is a calibrated parameter. This assumption is supported by empirical evidence (Duchemin et al., 2006;Groeneveld et al., 2007). We can therefore re-write the equation for 225 evapotranspiration in the form of a nonlinear regression: where ET i,s,t is average MODIS ET for an individual pixel i in season s and year y, T A,i,s is the interpolated post-2000 seasonal average air temperature (See Sect S3), and N DV I i,s,y is the average Landsat NDVI for the same year-season-pixel combination. The regression coefficients a s represent the estimated average effect of extra-terrestrial radiation R a and diurnal temperature range T R , and c mediates the effect of NDVI on the the crop coefficient. These coefficients are assumed to be 230 stationary across pixels (homogeneous) and years (stationary). We used a cross-validation analysis to evaluate how robust our ET estimates were to uncalibrated data. A calibration sample was formed by independently drawing 80% of pixels from each image, which we used to to estimate the regression coefficients a s and c. The estimated coefficients were then used to predict ET on the remaining 20% of the pixels using Equation 2. Predictions matched ET observation on the validation with a high degree of accuracy (R 2 = 0.87).

235
The cross-validation analysis results allowed us to build confidence in our assumptions that c and a c were stationary and homogeneous. We relied on these assumptions to predict pre-2000 ET using Equation 2 and regression coefficients a s and c estimated using post-2000 observations.
Finally, Strahler et al. (1999) provide annual classification of 17 land use categories using MODIS imagery. We combined these categories into six super classes to represent major land use categories within the watershed: grassland and shrubs, forest, The Upper Jhelum contains a variety of storage reservoirs including surface water bodies, groundwater, snow, and glaciers ( Figure 1 d-g). From the water balance, a decrease in storage (e.g., glacial melt) must be matched by a corresponding increase in the outgoing fluxes, evapotranspiration or streamflow. As such, a long-term decrease in catchment storage could have produced 250 a corresponding increase in streamflow followed by a similar decrease as catchment storage stabilizes. We hypothesize that there could have been a long-term reduction in permanent water storage leading to a subsequent decline in streamflow. For instance: -Hypothesis 6: A long-term decline in glaciers produced an increase and subsequent decrease in streamflow -Hypothesis 7: A long-term decline in permafrost produced an increase and subsequent decrease in streamflow 255 These hypotheses are grounded in studies of other catchments showing that a warming climate could temporarily increase streamflow through glacial loss (Singh and Kumar, 1997;Schaner et al., 2012) or permafrost melting (Kurylyk et al., 2016;Qiang et al., 2019). Such changes have been predicted to occur in high-elevation montane regions and could have contributed to the increase in streamflow in the 1980s and 1990s and subsequent decline after 2000 ( Figure 1c).
In addition to these long term effects, catchment storage at the seasonal time scale scale creates a time lag between precipi-260 tation and streamflow. Understanding how snow cover, groundwater, and surface water storage have changed over time would provide additional insight into the processes governing hydrological change in the catchment. Specifically, we hypothesize that: -Hypothesis 8: Reduced snow cover and earlier snow melt generated an earlier peak in annual hydrograph.
-Hypothesis 9: Reduced groundwater recharge led to a reduction in the baseflow contribution to streamflow throughout 265 the watershed.
All things being equal, earlier snow melt would tend to increase the amount of streamflow early in the year and decrease streamflow later in the year. Earlier snow melt would also allow greater vegetation activity and evapotranspiration earlier in the spring, thereby reducing groundwater recharge. In aggregate, the processes involved in these two hypotheses might have increased the proportion of annual precipitation exiting the catchment as evapotranspiration instead of streamflow. We also consider the potential for loss of permafrost to have produced an increase and subsequent decrease in streamflow.
For example, Qiang et al. (2019) found that melting permafrost generated a temporary increase in streamflow in the upper Yellow river of 5%. To evaluate the possibility for this process in the Upper Jhelum, data were downloaded from the Global Permafront Zonation Index (GPZI) map (Gruber, 2012). Given the uncertainty in permafrost occurrence, the GPZI is presented 280 on a scale that indicates the likelihood of permafrost, with a minimum indicating that "permafrost exists only in most favorable conditions" and maximum indicating that "permafrost exists in nearly all conditions." We binned this scale into five groups of permafrost likelihood including low, medium-low, medium, medium-high, and high. The upper Jhelum contains no pixels with medium-high or high likelihood of permafrost and in most of the areas where permafrost is possible, the likelihood is low (see Fig. S8). To evaluate the potential for permafrost to affect streamflow, we compared the areal extent of the GPZI with the 285 necessary loss of frozen water storage to produce the observed changes in streamflow.

Data sources: Snow
Winter precipitation occurs largely as snowfall and remains in some parts of the catchment until the late summer. Because different regions of the watershed may be affected by missing pixels (e.g., clouds) on any given acquisition date, we separated the watershed into 15 distinct zones of roughly equal areas defined by three elevation bands (<1650 m, 1650-2200 m, and 290 >2200 m) in the five local subwatersheds corresponding to the available stream gauges. Snow contains a distinct spectral signature with high reflectance in visible and near-infrared bands and low reflectance in shortwave infrared bands, and can therefore be detected by from normalized different snow index (NDSI), which is defined as (Green -SWIR)/(Green + SWIR) (Dietz et al., 2012). We generated timeseries of snow cover in each of the 15 zones using Landsat 5 imagery by applying a threshold of 0.5 to NDSI to distinguish snow and water cover from dry land. We further distinguish snow (bright) from open 295 water (dark) using a threshold of 0.2 on the NIR band reflectance (Kulkarni et al., 2002). For each zone, we selected only the dates where missing pixels constituted less than 25% of the zone, leaving an average of 79 and 65 observations in each zone before and after 2000, respectively. We generated an adjusted snow estimate for each zone and date by dividing raw snow cover estimates by the fraction of missing pixels. We then estimated average daily snow cover by fitting locally weighted non-parametric time regressions (LOESS) (Cleveland et al., 1992) to snow cover observations, for each zone, before and after 300 2000. The final smoothed estimate of snow cover was taken as the sum of snow cover across all zones, with 95% confidence intervals calculated from the sum of standard errors across the zones.

Data sources: Groundwater
We were unable to obtain in situ groundwater observations and remotely sensed observations from GRACE satellite were inadequate due the large spatial averaging kernel (≈ 40,000 km 2 , compared to a catchment area of approximately 13,000 km 2 ) 305 and lack of observations prior to 2002. We therefore relied on proxies to assess changes in groundwater behavior.
First, we considered changes in baseflow discharge an indication of changes in groundwater. Commonly applied (non-)linear reservoir models assume that baseflow discharge follows a monotonic relationship with groundwater storage. Therefore a change in baseflow over time (see Sect. 2.2) would indicate a change in groundwater storage.
Second, we looked for hysteresis in the relationship between surface water storage and streamflow as indicative of a ground-

Data sources: Surface water
A number of lakes and wetlands exist throughout the valley including Wular lake, which intersects the main stem of the Upper Jhelum between gauges (iii) and (iv), and seasonally inundated valley wetlands which capture flow from the subwatershed that 320 drains into gauge (iv) (see Fig. 1). The actual volumetric surface water storage of the catchment is difficult to estimate. Instead, we focus on changes in surface water area using remote sensing imagery. We classify surface water extent in Wular lake and in the wetlands in all available Landsat imagery over the period 1984-2013.
Open water is highly absorptive in short-wave infrared bands and more reflective in bands with shorter wavelengths. We use the modified normalized-difference water index (MNDWI Xu, 2006) as an indication of the likelihood of open surface 325 water, with a threshold distinguishing between land and water pixels. Because water exhibits spatial coherence dictated by topography, we gap-filled missing pixels in Landsat 7 bands due to the scan-line corrector error by propagating edge pixels towards the center of the gap (as in Penny et al., 2018). We used a fixed MNDWI threshold across all images to classify surface water. Clouds were identified using the rudimentary Simple Cloud Score algorithm in Google Earth Engine (ee.Algorithms.Landsat.simpleCloudScore()) on top-of-atmosphere images. We applied this classification ap-

Results
Our estimates of rainfall, streamflow, evapotranspiration and storage change allowed the water balance to be closed for the 1984-1999 and 2000-2013 periods with an error of ±15% (Fig. 3). This suggests that the total watershed fluxes are estimated with reasonable confidence, particularly given uncertainty in rainfall interpolation and remote sensing models of evapotranspi-335 ration.
However, the analysis did not allow us to close the differential water balance (i.e., changes in the water balance) between the two periods. We observed an average decrease in precipitation of 117 mm per year and an increase in ET of 32 mm per year, to attribution does not rely on our ability to close the differential water balance. Instead, the hypotheses developed above (Section 3) are investigated individually, as described below. While data limitations might prevent a subset of the hypotheses from being conclusively tested, the approach nevertheless reveals a process-based understanding of hydrological change in the Upper Jhelum watershed.

Precipitation: hypotheses 1-3 350
Precipitation exhibited notable changes in total volume, as annual precipitation fell by 117 mm, corroborating Hypothesis 1.
Precipitation also exhibited changes across seasons (Fig. 4), consistent with Hypothesis 2. These changes were driven almost entirely by a loss of spring precipitation of 117 mm. The other seasons saw modest and statistically insignificant changes in precipitation (+20 mm in winter, -14 mm in summer, -5 mm in autumn).
Additionally, the number of storms greater than 1 mm increased in all seasons except spring, thus supporting Hypothesis 3 355 (Fig. 4c). Combined with the fact that precipitation either decreased or remained constant, this entails a decrease in average Overall, of the 32 mm annual increase in ET that we detected, approximately 17% can be attributed directly to increasing air temperature, and the remaining 83% to an increase in NDVI. The largest net contributors to watershed-averaged increases in 380 ET were shrubs and grassland (15 mm), cropland (11 mm) and forest (3 mm Overall, roughly 0.19% of the watershed has at least a medium likelihood of permafrost (Fig. S8). In a larger portion (6.4%) permafrost has a low likelihood and therefore requires extremely favorable conditions and is therefore likely to be shallow (Gruber, 2012). In that context, the 50 mm per year increase in Upper Jhelum streamflow necessary to close the differential water balance would require a loss of 780 m per year of permafrost in all regions where permafrost is likely to be shallow. Such Snow cover follows a seasonal signal, with a maximum in late winter followed by receding snow pack in spring. After 2000, the maximum snow cover cover extent was reduced and was followed by more rapid snow melt (Fig. 7a). The extent of 405 the influence of snow on the hydrograph appears to be reflected by an earlier peak in the hydrograph. This finding supports Hypothesis 8, although observational data is missing to quantify the implicated changes in storage volumes. Earlier snow-melt and its hypothesized effect on the water balance is nonetheless consistent with the observed increased in ET during the spring (see Sect. 5 for further discussion).
Remote sensing observations of open water extent suggests that surface water storage declined dramatically during the study 410 period, both in Wular lake and the neighboring wetlands (Fig. 7b). Both of these surface reservoirs connect to the main stem of the Upper Jhelum between gauges iii. and iv. (see Fig. 1), and likely play an important role in streamflow generation along this reach. To investigate this effect, we computed locally generated baseflow as the difference in baseflow between gauges iii. and iv. Baseflow along this reach peaks earlier and at a much lower amplitude after 2000 (Fig. 7c). It also transitions much earlier (mid-summer) to losing conditions (i.e. negative local streamflow values), compared to pre-2000 where baseflow transitioned to 415 losing conditions only at the end of autumn and beginning of winter. Additionally, the relationships between locally generated baseflow and open water extent (of the lake and the wetlands) exhibit a seasonal hysteresis (see Fig. S7), which we interpret as indication that groundwater plays an important role mediating streamflow generation in the valley bottom. Although the shape of the hysteresis remains similar, the amplitude of the hysteresis behavior is reduced considerably after 2000.
The implication of changing groundwater storage can be better understood by comparing temporal trends in streamflow 420 below Wular lake (iv., Sopore station) and the most upstream gauge of the watershed (i., Sangam station). The streamflow timeseries at both locations exhibit statistically significant decreasing trends over the period 1960-2013 (Fig. 8a). Baseflow separation, however, reveals important differences between both streamflow timeseries. Baseflow decreases over time only in the downstream gauge (Fig. 8b) whereas quickflow decreases only in the upstream gauge (Fig. 8c). Temporal changes in the baseflow index (B = Q B /Q T otal ) of each of these gauges therefore occurs in opposite directions, with decreasing B near the 425 outlet and increasing B in the hinterlands (Fig. 8d). This lends credence to Hypothesis 9, and we further discuss potential causes and implications of these opposing trends in the baseflow index with respect to saturated and unsaturated groundwater storage in the Discussion (Sect. 5).

Discussion
We now synthesize the results presented above to develop a narrative of hydrological change by reconciling the various fluxes 430 that have either increased or decreased over time (Fig. 9, red and blue arrows). Taken together, the observed changes in hydrological fluxes indicate additional unobserved changes in fluxes connecting surface and groundwater that might play an important role in explaining hydrological change in the catchment. As discussed in the following paragraphs, evidence suggests that these fluxes have decreased over time (Fig. 9, pink arrows).
-The largest seasonal precipitation occurs in spring, when the prevailing climate is driven by westerlies, yet spring pre-435 cipitation declined considerably during the study period. At the same time, vegetation activity increased across most of the watershed within both anthropogenic (cropland, orchards) and natural (forest, shrubs and grassland) land use classes (Fig 5a), producing an increase in evapotranspiration. The corresponding reduction in spring streamflow was partially compensated by higher temperatures and earlier snow melt. Near the end of the study period, the snow pack was nearly always exhausted by the end of spring (Fig 7a). -Springtime seepage generates high groundwater recharge throughout the valley, both from the highlands and from within the valley (Jeelani, 2008). High groundwater levels at the end of the spring season are reflected by the peak in baseflow along the reach of the stream that includes Wular lake (Fig. 7d). Prior to 2000, this peak occurred near the transition between spring and summer. After 2000, this peak occurred earlier in spring and at a much lower level of baseflow. This suggests considerably lower groundwater storage over time, reflected in the reduced groundwater recharge fluxes (Fig 9,   445 Spring).
-In summer, the prevailing climate is driven by the Indian monsoon and precipitation is generally less than westerly precipitation in the spring (Fig 4a). Consequently, hydrology within the watershed is controlled largely by water storage (snow, lakes and groundwater) left over from spring. Prior to 2000, the seasonal recession of the baseflow hydrograph starts high at the end of spring and continues to produce discharge throughout the summer season (Fig 7d). The river 450 then transitions to losing conditions in autumn and earlier winter before baseflow increases with winter precipitation and snow melt (Fig 7d). After 2000, snow pack at the beginning of summer is nearly exhausted, and surface water (Fig. S5) and groundwater (Fig. 7d) storage are also reduced. The receding limb of the hydrograph starts low at the beginning of summer and quickly depletes, transitioning to losing conditions within the summer season. Indeed, the recovery of local baseflow in autumn (i.e., baseflow becoming less negative, see Fig. 7d, brown) is driven by declining baseflow 455 downstream rather than increasing baseflow upstream. Additionally, summer evapotranspiration after 2000 increased throughout much of the watershed including in natural and anthropogenic land use classes.
-The hydrological cycle in the watershed is mostly dormant in autumn and winter. In autumn, little precipitation falls and the hydrograph recedes into winter. Notably, lake storage depletes and the lake transitions to losing conditions in several years, particularly before 2000 (Fig. 7d). Winter precipitation arrives mostly as snowfall, replenishing snow storage.

460
Winter rain and snow melt serve as the early primers for the seasonal cycle to renew in spring.
To summarize, climate appears to be the primary cause of hydrological change within the Upper Jhelum. The most influential driver is the decline of spring westerly precipitation. Other studies have associated this decrease with warming temperatures (Zaz et al., 2019). This effect is compounded by an array of other drivers that affect watershed processes. Notably, the loss of baseflow downstream of Wular lake suggests a decrease in groundwater storage in the valley. This decline in groundwa-465 ter is facilitated both by reduced rainfall during the spring and increasing watershed evapotranspiration. The latter might be exacerbated by an increase in the number of storms for precipitation falling outside of the spring season. This observation from the PERSIANN CDR dataset allows us to hypothesize that the shift towards a larger number of smaller storms results in reduced overland and macropore flow, along with more stable and persistent soil moisture and ultimately more water "lost" to vegetation uptake. In this case, seepage would be increasingly likely to occur via slow drainage processes, rather than macro-470 pore flow activated in large storms. Such changes have been observed in other karst catchments (Zhao et al., 2019) and are supported by the evidence that quickflow declined and the baseflow index increased in the most upstream gauge in the Upper Jhelum (Fig. 8). The increase in evapotranspiration appears to have occured throughout the catchment. Our evapotranspiration model indicates that increasing air temperature had a small direct effect on ET (a 2% increase) and that most of the overall The climate signal therefore manifested itself by providing favorable conditions for plant biomass growth, as represented by increasing greenness. Such a change would be expected with warming temperatures and reduced snow cover early in the year, in addition to greater energy availability (temperature) and water availability (soil moisture) that arose due to climate warming and increased storm frequency, respectively. Evapotranspiration exhibited the greatest increases in regions that transitioned to orchard cultivation, but these areas represent a small fraction of the overall watershed and increase in ET. Evapotranspiration 480 also increased considerably in cropland areas, likely due to the aforementioned climatic changes in addition to agricultural intensification. Indeed, increasing fertilizer use in the catchment has been noted in other studies (Wetlands International South Asia, 2007) and might have contributed to the increased ET in agricultural land use classes. It is also an important driver contributing to the growth of aquatic vegetation and increasing ET in Wular lake in the summer (see Figs 5,6).

485
In this study, we develop a bottom-up approach to hydrological attribution to understand the drivers of dramatic changes in the annual streamflow of the Upper Jhelum river. We found that much of the observed decrease in streamflow is associated with decreases in westerly precipitation in spring, in addition to greater evapotranspiration. While land-use change to orchard plantations and agricultural intensification are likely contributing factors, we attribute most of the increase in evapotranspiration to non-local anthropogenic causes, most notably increased vegetation activity in spring due to increased temperature and earlier 490 snow melt.
The approach focuses on the development of a process understanding of hydrological change in the catchment by separately evaluating multiple hypotheses about specific mechanisms, rather than relying on a basin-wide aggregate model for predictive inference. As such, it is a promising tool to attribute hydrological change in situations where process uncertainty might be compounded by hydrologic regime shifts. Additional uncertainty in the water balance, due to bias in hydrological observa-495 tions, limits the utility of catchment-scale models, including a model as simple as the water balance. Indeed, detection biases associated the different remote sensing sources might compound and prevent the proper calibration of an aggregate model, as seen in our inability to close the differential water balance. Instead of being used as a predictive tool, the water balance (or any other conceptual model) is here used as a framework to generate hypotheses that can be investigated individually. While outcomes associated with individual hypotheses might exhibit considerable uncertainty, especially in data-scarce catchments, 500 together the multiple hypotheses provide multiple strands of evidence to support (or refute) specific mechanisms and ultimately attribute hydrological change.
Data availability. The secondary data that support the findings of this study are available on request from the corresponding monitoring and collection agencies. Remote sensing products used in this study are freely available from the respective providers (see in-text citations).