Technical note: Accounting for snow in the estimation of root-zone water storage capacity from precipitation and evapotranspiration fluxes

A common parameter in hydrological modeling frameworks is root-zone water storage capacity (SR[L]), which mediates plant-water availability during dry periods and the partitioning of rainfall between runoff and evapotranspiration. Recently, a simple flux-tracking based approach was introduced to estimate the value of SR (Wang-Erlandsson et al., 2016). Here, we build upon this original method, which we argue may overestimate SR in snow-dominated catchments due to snow melt and evaporation processes. This matters for hydrological predictions under warming scenarios with decreased snowpack. 5 We propose a simple extension to the method presented by Wang-Erlandsson et al. (2016), and show that the approach provides a more conservative minimum estimate of SR in snow-dominated watersheds. This SR dataset is available at 1 km resolution for the continental United States, along with the full analysis code, on Google Colaboratory and Earth Engine platforms. We highlight differences between the original and new methods across the rain-snow transition in the Southern Sierra Nevada, California, USA. As climate warms and precipitation increasingly arrives as rain instead of snow, the subsurface may be an 10 increasingly important reservoir for storing plant-available water between wet and dry seasons; improved estimates of SR will therefore better clarify the future role of the subsurface as a storage reservoir that can sustain forests during seasonal dry periods and episodic drought.

plants: there is considerable evidence of plant water uptake from depths below soil (Dawson et al., 2020;Schwinning, 2010).
For example, roots can extend into and draw water from the bedrock vadose zone (rock moisture, sensu Rempe and Dietrich, 2018;Hahm et al., 2020) or groundwater (Miller et al., 2010;Lewis and Burgy, 1964). Within seasonally dry environments in 25 particular, a significant volume of water accessed during the growing season can be derived from depths below mapped soils (Rose et al., 2003;Jones and Graham, 1993;Arkley, 1981). We emphasize that an accurate representation of S R therefore should include not only moisture available within the soil, but also plant-accessible water below the soil.
S R does not, however, include snowpack, which is an above-ground water storage reservoir. Conservatively estimating S R in systems that currently receive a significant proportion of their precipitation as snow is particularly important given the 30 ongoing shift from snow to rain under a warming climate, and the attendant heightened significance of subsurface water storage dynamics to plant ecosystems and streams. An existing, widely used method for estimating S R (Wang-Erlandsson et al., 2016) does not account for snow melt and sublimation/evaporation dynamics, which may result in significant overestimation of this reservoir. We present an extension to the original method to ensure a conservative estimate of S R in snowy areas. We describe the method details and highlight results from a rain-snow transition transect in the Southern Sierra Nevada, California, USA. 35 We also provide a geotiff raster map of S R across the continental United States at the 1 km pixel scale. Finally, we link to a Google Earth Engine (https://earthengine.google.com/) script written in Python (https://www.python.org/) within the Colab coding environment (https://colab.research.google.com/) to document application of the method, and to facilitate comparative analyses using other widely available and spatially distributed precipitation, snowcover, and actual evapotranspiration datasets.

40
To estimate S R , Wang-Erlandsson et al. (2016) compute a running root-zone storage deficit (more positive means larger capacity in the subsurface for storage) using differences between fluxes exiting (F out ) and entering (F in ) the root zone during a given time interval (typically equal to the sampling period of the remotely sensed evapotranspiration dataset). To obtain a conservative estimate of S R (that is, a robust lower bound), it is important to make sure that F in is not underestimated (when in doubt, assume all precipitation enters the rooting zone), and that F out is not overestimated (when in doubt on the amount 45 of F out that contributes to increases in the root zone storage deficit, simply set F out = 0). Accordingly, other (often difficult to measure at large spatial scales) fluxes exiting the rooting zone, such as leakage and runoff, are ignored when computing running storage deficits, as neglecting these fluxes results in a conservative estimate of S R . Due to these simplifications, F in and F out are set equal to precipitation (P ) and evapotranspiration (ET ), respectively.
The original storage deficit tracking (and subsequent estimation of S R ) procedure presented by Wang-Erlandsson et al.

50
(2016) is achieved through two steps. First, over a given time interval t n to t n+1 , the accumulated difference (A tn→tn+1 ) between F out and F in is calculated as: Here, since the root-zone storage deficit is being calculated (and not actual storage), the incoming and outgoing fluxes have opposite signs from a conventional mass balance. A lower bound on the root-zone storage deficit at each time interval can then 55 be calculated as the maximum value of 0 and the running sum of these accumulated differences: Finally, S R is estimated as the maximum observed value of D.
The potential inaccuracy introduced by this original method that we explore here is that during periods when snowpack is present within the pixel, F in may be underestimated (due to melting snow entering the rooting zone, for example, when 60 precipitation is zero) and F out from the root zone may be overestimated (due to attribution of sublimation/evaporation from the snow surface to a flux from the subsurface). As discussed above, both of these possibilities may lead to overestimation of In the absence of spatially and temporally resolved information about snowmelt and sublimation dynamics, a simple and conservative way to correct for these potential errors is to continue to decrease the storage deficit as incoming precipitation 65 arrives, and to set F out = 0 during periods when snow cover (C, the fraction of the pixel covered in snow, which is reliably measured at large spatial scales via satellites) is present, thereby not counting evapotranspiration towards increasing the storage deficit during snowy periods. We therefore introduce a correction term for the outgoing flux in the calculation of the accumulated difference between outgoing and incoming fluxes during each time interval: where C 0 is some threshold below which it is assumed that snow cover is negligible, and · is the ceiling operator (rounding up to the nearest integer), returning a 1 if C > C 0 and 0 if C ≤ C 0 . The expression therefore effectively sets F out = 0 whenever snow is present (or deemed not-negligible) in the pixel, providing a robust lower-bound estimate of S R in the running storage deficit calculation.

75
We implement the original and snow-corrected algorithm developed here using Google Earth Engine, accessed via a Google Colab notebook and the Python programming language's Earth Engine application programming interface. This readily enables i) access to distributed timeseries hydrological products (i.e., snowcover, evapotranspiration, and precipitation); ii) computation in the cloud and iii) a shareable script that can be quickly modified and executed by new users (see link at end of manuscript).
The algorithm requires precipitation and evapotransipration datasets to compute F in and F out , and a snow cover dataset to 80 implement the proposed snow-correction step. We use Oregon State's PRISM daily precipitation product (https://developers. , and set C 0 = 0.1 (snow cover is assumed negligible at less than C 0 = 10% pixel coverage; in this case, C 0 = 10% is also the resolution of the underlying snow cover dataset). We restrict our analysis to the temporal intersection of

Results
To demonstrate the impact of the method in snowy areas, we focus on a region of the Sierra Nevada in California, United States, where there exists a strong elevation-driven gradient in winter snow cover. Figure 1 illustrates three raster data layers derived from application of the new method. Figure 1a plots root-zone storage capacity calculated using the snow-correction method.

105
Values range from near 0 mm over exposed bedrock outcrops in the High Sierra, to over 900 mm in the dense mid-elevation forests. Figure 1b shows the difference between S R computed using the original method and the snow-corrected S R . Figure   1c plots average winter (January through April) snow cover. As expected, the difference in Figure 1b is small in the lower, rain-dominated elevations, and larger in areas with snow cover. However, some areas with substantial snow cover show small differences between the methods. These are likely areas where S max is small, coinciding with low-energy, exposed-bedrock 110 locations at high elevations. Figure 2 illustrates the full timeseries output of the snow-accounting and original methods at two locations, identified by white points in Figure 1. The location farther west is a 'low snow' location, with negligible snowfall (snow present less than 1% of the time) during the winter months, and the location to the east is a 'high snow' location, with snowcover present over 50% of the time during the winter months. Gray shading in all subplots indicates that greater than 10% of the pixel is covered  in snow at that time point, during which evapotranspiration is set to zero in our method (lower panels). The top panels of Figure 2 plots storage deficits using the original and the snow-accounting methods, clearly demonstrating the divergence of deficit calculations between the two methods in the region with significant snow cover. In all instances, S R is calculated as the maximum observed value of the storage deficit. In the high snow location using the original method, this leads to an estimated value of S R that is approximately 50% larger than that calculated with the snow-accounting method.

Discussion
Our proposed method for estimating S R is a minimum estimate. Actual S R should generally exceed estimated S R values presented in our revised method, because some evapotranspiration occurs during times when snow cover is present. The snowaccounting method and the original method do not account for leakage, surface runoff, and upslope drainage in the calculation of F in .

125
S R in rain-dominated climates has been shown to impact drought resilience (Hahm et al., 2019a), and snow-rain transition elevations are increasing as the climate warms (Knowles et al., 2006). If precipitation arrives as rain rather than snow, the role of the subsurface in storing that water for plants will likely be amplified. Globally, mountainous snow-rain transition zones tend to coincide with forested areas, underscoring the importance of accurate estimates of S R for prediction of forest sensitivity to climate variability in the future.

Conclusions
We argue that an existing method for estimating root-zone water storage capacity (S R ) will tend to overestimate S R in snowy areas due to unaccounted for snow melt, evaporation, and sublimation processes. We provide a correction factor that relies on a widely available distributed percent snow cover dataset to provide a tighter lower bound estimate on S R . Conservatively