Detecting Streamflow in Dryland Rivers Using CubeSats.

Determining the flow regime of non-perennial rivers is critical in hydrology. In this study, we developed a new approach using CubeSat imagery to detect streamflow presence using differences in surface reflectance for areas within and outside of a river reach. We calibrated the approach with streamflow records in the Hassayampa River of Arizona over 3 years (2019-2021), finding good agreement in the annual fractions of flowing days at stream gages (R 2 = 0.82, p < 0.0001). Subsequently, annual fractions of flowing days were derived at 90 m intervals along the Hassayampa River, finding that 12% of reaches were classified as intermittent, with the remaining as ephemeral. Using a Hovmöller diagram, streamflow presence was visualized in unprecedented spatiotemporal detail, allowing estimates of daily fraction of flowing channel and annual fractions of flowing days. This new tool opens avenues for detecting streamflow and studying hydrological and biogeochemical processes dependent on water presence in drylands.

Introduction

Intermittent rivers and ephemeral streams (IRES), which do not have continuous water throughout the year, are globally prevalent (Busch et al., 2020; Messager et al., 2021), especially in arid and semiarid regions, often referred to as drylands (Hammond et al., 2021). IRES within global drylands play critical roles as sources of water supply (Levick et al., 2008; Tooth, 2000), as habitats for plant and animal species (Datry et al., 2014; De Jong et al., 2015; Stromberg & Merritt, 2016), as hot spots for biogeochemical and contaminant reactions (Brooks et al., 2006; McKenna & Sala, 2018; Reneau et al., 2004), and as corridors for the sediment transport (Sims, 2010), surface‐groundwater exchanges (Vivoni et al., 2006), and stormwater effluent (Bateman et al., 2015). Naturally, the streamflow presence in IRES occurs under high degrees of spatial and temporal variability (Goodrich et al., 2018; Shanafield et al., 2021). Under the combined effects of climate change and human activities, the spatial and temporal extents of IRES have been increasing in recent years (Chalise et al., 2021; Price et al., 2021; Sauquet et al., 2021), and the trend is expected to continue (Döll & Schmied, 2012; Jaeger et al., 2014).

Despite growing scientific evidence, policymakers often do not recognize the importance of IRES in providing ecological and hydrological functions in landscapes (Acuña et al., 2014; Doyle & Bernhardt, 2011), and these river reaches usually lack effective protections (Datry et al., 2017; Skoulikidis et al., 2017). In 2020, for instance, changes to the definition of Waters of the United States removed the protection of ephemeral rivers from the Clean Water Act and raised concerns on water quality degradations (Fesenmyer et al., 2021; Mazeika et al., 2019; Mihelcic & Rains, 2020). Changing policies also challenge environmental agencies to determine the jurisdictional status of IRES for surface water protection (Taylor et al., 2011). Current efforts often include the use of streamflow gaging stations, field inspections, and citizen science reports which have low sampling capacities and are difficult to scale up to large regions.

A promising new technology is offered by small satellite remote sensing which has been used for surface water detection (McCabe, Aragon, et al., 2017, McCabe, Matthew, et al., 2017). For instance, Planet Labs uses a constellation of CubeSats (i.e., a large fleet of small satellite units) to survey the entire Earth at a high spatial resolution (∼3 m) and an unprecedented daily frequency since 2017. The resulting PlanetScope imagery have been used in a number of scientific applications in the natural and built environments (Houborg & McCabe, 2016; Kindler et al., 2022; Shugar et al., 2021; Vivoni et al., 2020), including for the detection of permanent and ephemeral water bodies and their changes over time (Cooley et al., 2017, 2019; Mishra et al., 2020; Panda et al., 2020).

Here, we developed a new approach for detecting rapidly changing streamflow in dryland rivers using PlanetScope imagery. Our robust method relies on surface reflectance contrasts of flowing water when present in arid and semiarid rivers as compared to their dry surroundings. Calibrated with streamflow gaging data available over short periods, the approach provides unparalleled spatial detail at near daily resolution. Its broad applicability is shown in a large region and a method is used to link its single parameter to flowing conditions. Its coverage over large remote areas offers scientists a new means to characterize non‐perennial rivers and environmental management agencies a novel tool in the jurisdictional determination of IRES.

Materials and Methods

Study Area

We apply the approach to the Hassayampa River (HR, 186 km in length) in Arizona (Figure 1a). The north–south oriented watershed covers ∼3,766 km2 of the Sonoran Desert. The elevation relief (240–2,430 m) leads to large changes in climate, vegetation, and flow regimes in the watershed (Mascaro, 2017; Medgal et al., 2008; Stromberg et al., 2009). Mean annual precipitation varies from <200 to 437 mm, with maxima in the summer (50%–58%, July to September) and winter (30%–39%, November to March) according to climate normals from 1991 to 2020. Land cover types, shown in Figure 1b, are shrublands (85.6%), forests (6.9%), grasslands (5.3%), urban areas (1.5%), and croplands (0.3%) in the National Land Cover Dataset (NLCD) 2019 product (Jin et al., 2019), and organized into two ecoregions: AZ/NM Mountains (ANM) and Sonoran Basin and Range (SBR, https://www.epa.gov/eco-research/ecoregions). ANM has higher elevations and receives more precipitation that feeds the perennial headwaters of HR (Ruhlman et al., 2012), while SBR is characterized by a warm and arid climate that is dominated by ephemeral streams (Calzia & Wilson, 2012). For jurisdictional evaluations of surface waters, the Arizona Department of Environmental Quality (ADEQ) divided the main stem of the HR at major confluences (Figure 1b). For each segment, ADEQ either assigned a flow regime (perennial, intermittent, or ephemeral), or classified as undetermined if there is no available data (http://azdeq.gov/FlowRegimes). For the nine subreaches along the HR, four have undetermined flow regimes, and the rest are classified as perennial or intermittent (Table S1 in Supporting Information S1).

Figure 1

Study area and its characteristics. (a) Location of the Hassayampa River in central Arizona. (b) Land cover classification from National Land Cover Dataset (NLCD) at 30 m resolution, with study reaches along the HR and streamflow gages operated by Flood Control District of Maricopa country (FCDMC) and US Geological Survey (USGS). Names of available gages are also labeled. The footprints of PlanetScope tiles are shown in blue boxes, with the tile IDs labeled for reference. The dashed line shows the boundary of two ecoregions: AZ/NM Mountains (ANM) and Sonoran Basin and Range (SBR). (c) Channel width arranged from upstream to downstream (gray lines represent raw data and black lines represent smoothed results), with the extent of reaches marked with different colored backgrounds and the ecoregion boundary (dashed line).

Ground Observations

The HR is monitored by a streamflow gage and weather station network from the Flood Control District of Maricopa County (FCDMC) and the US Geological Survey (USGS), shown in Table S2 in Supporting Information S1. For each streamflow gage, we obtained precipitation records from a collocated rain gage (available at WAG, BCN, and U60 sites) or the nearest rain gage (<5 km in distance), all from FCDMC. The paired precipitation and streamflow records at six sites for the water years (WY) of 2019–2021 (October 2018 to September 2021) are described in Text S1 in Supporting Information S1 to provide a first indication of spatiotemporal variations in precipitation and streamflow along the HR.

PlanetScope Imagery From CubeSats

PlanetScope data are acquired by a constellation of CubeSats operated by Planet Labs. The system provides data in the visible (Red, Blue, and Green) and near infrared (NIR) bands with 3‐m resolution (Planet Labs, 2021a). Here, we used merged, radiometrically‐, sensor‐, and geometrically‐corrected PlanetScope imagery with defined grids organized into 25 km by 25 km tiles covering the main stem of the HR (labeled as 909, 809, 808, 708, 608, 508, and 408 in Figure 1b) of three WYs (2019–2021). All images were surface reflectance products, which were atmospherically corrected using the 6SV2.1 radiative transfer model (Planet Labs, 2022). Details on the PlanetScope data processing and availability are described in Text S2 in Supporting Information S1.

Streamflow Detection Method

The streamflow detection method is based on robust differences in surface reflectance for areas inside and outside of a river reach (see workflow in Figure S3 in Supporting Information S1). To illustrate the method, Figure 2 shows an example at the U60 site for conditions prior to and during a flood in the HR. We derived the average surface reflectance values for pixels inside and outside the river reach containing the gage site (see Text S3 in Supporting Information S1 for details on buffer zones and the reach mask). We tested various bands and band ratios (Figure S4 in Supporting Information S1) and identified that the NIR band exhibited the largest differences between flowing water inside a reach and dry conditions in the surrounding buffer zone. This selection is supported by the strong absorption properties of water in the NIR band (Mishra et al., 2020). For instance, NIR in the reach (NIRIn) significantly decreased from dry to flood conditions (Figures 2d and 2e), while NIR outside of the reach (NIROut) showed much lower changes. Based on this, we used the difference between NIRIn and NIROut in each buffer zone or as an index to determine flowing conditions. Since each buffer zone is 90 m in river length, a high spatial resolution of NIRDiff can be achieved (Figure 2f) and compared with streamflow gaging data. We determined a threshold value of NIRDiff that best matched the fraction of flowing days (FF ) obtained from the approach (F NIR) and from stream gages (F Gage) over coincident times when usable PlanetScope imagery were available (in the following, all fractions are shown as percent values). In this study, a NIRDiff threshold (θ) of −0.05 was selected using all gages (see detailed descriptions in Text S4 in Supporting Information S1). Given the representativeness of the gaging sites, this value of θ was applied along the entire HR. To determine the flow regime (i.e., whether a reach is ephemeral or intermittent) along the entire HR, we compared the average annual FF from the approach to previously‐determined values for different ecoregions by converting flowing days thresholds from Granato (2010) to FF values. For the two ecoregions in the HR, a river reach with FF  ≥ 7.1% in the ANM ecoregion (FF  ≥ 3.6% for SBR) is considered as intermittent. Otherwise, the flow regime is classified as ephemeral.

Figure 2

Remote sensing method and its calibration with field data. (a–b) True color images from PlanetScope under dry (03/09/2020) and flood (03/20/2020) conditions at the U60 site (image is 500 m in radius centered on gaging station shown as a triangle). (c) Reach mask and buffer zones at the U60 site. (d–e) Surface reflectance of NIR under dry and flood conditions. (f) ΔNIR computed as flood minus dry condition. (g) Scatter plot of the annual fraction of flowing days from gages (F Gage) against the annual fraction of days with ΔNIR lower than the optimal threshold (F NIR). Different colors in the legend correspond to different gages. Dashed line is 1:1 line and solid line is the linear regression, with additional details presented in Text S4 in Supporting Information S1.

Results

Temporal Dynamics at a Gaging Station

We tracked the temporal dynamics of river conditions at the HR gage stations using ground and PlanetScope data for WYs 2019 to 2021. For instance, Figure 3a presents daily precipitation (P) and streamflow discharge (Q) at the U60 site. During the 3‐year study period, annual total P was slightly below average (95% of 1995–2020 average at the site), but the seasonal distribution exhibited lower summer P and higher winter P than average (65% and 111%, respectively). Over this period, PlanetScope data captured NIR differences between locations within and outside the channel reach (Figure 3b). When the channel was dry, for example, from May 2020 to June 2021 with no streamflow, NIRIn was consistently higher than NIROut (0.34 and 0.27, respectively). This was due to the higher reflectivity of the dry substrates in the active channel as compared to those outside and adjacent to the reach (NIRDiff of +0.07, Figure 3c). Similar patterns also occurred during other dry periods (e.g., April–August 2019). When flowing water was present, however, NIRIn crossed below the relatively stable NIROut since the surrounding areas had limited changes in NIR, resulting in a negative NIRDiff. Drying of the reach led to increases in NIRIn which eventually became larger than NIROut, creating a second cross‐over point. The duration of the cross‐over period with NIRDiff threshold θ ≤ −0.05 matched the fractions of flowing days (FF ) obtained at the gaging station (blue shaded periods in Figure 3c), such that the method is a robust way of determining streamflow presence.

Figure 3

Temporal dynamics at a streamflow gaging station. (a) Daily total precipitation (P) and daily average streamflow discharge (Q) at the U60 site during the study period (WY 2019–2021). (b) Raw (dots) and smoothed (lines) NIR inside the reach and outside the reach mask within the 90 m buffer segment containing the U60 site. (c) Raw (dots) and smoothed (line) NIRDiff as compared to the flowing days from the gaging records (shaded blue areas). Streamflow detection occurs when NIRDiff is below the threshold (θ) of −0.05 (dashed horizontal line).

Streamflow Presence Along the Hassayampa River

We then applied the method to the entire main stem of the HR and summarized the results using a Hovmöller diagram of NIRDiff organized by river distance and date (Figure 4a). Each pixel in the diagram represents the average NIRDiff of each 90 m buffer zone for each day. Lower NIRDiff (i.e., blue colors) indicate streamflow presence, while NIRDiff values near or above zero (i.e., yellow colors) show periods and locations with no water. Pixels with no data are shown in gray, either due to the lack of high‐quality imagery (vertical features) or the inability to determine the reach (horizontal features). The Hovmöller diagram allowed a rapid assessment of the spatiotemporal presence of water along the HR during the entire record. Clearly, the upper HR reaches (Reach 7B, 5, 4, and 3) had more flowing days than reaches at lower elevations (Reach 2B and 1A), except for Reach 1B where agricultural return flows led to frequent water presence. Periods of time with streamflow presence in the winter and summer seasons were also evident from the diagram. This was summarized using the daily fraction of flowing channel (FFC, Figure 4b) calculated as the fraction of the HR with streamflow presence in the reach mask. Low FF (average of 2.5% and standard deviation of 3.5%) indicated that HR was non‐perennial overall but had strong seasonal variations. FF values of winters were usually higher than summers, except for WY 2021 when summer storms led to streamflow presence in many reaches (e.g., FF  = 31% and 33% on July 18 and 19 August 2021).

Figure 4

Spatiotemporal dynamics of streamflow presence. (a) Hovmöller diagram of NIRDiff as a function of distance and date. Gray areas depict missing data due to no usable PlanetScope imagery (vertical features) or undetectable reaches due to vegetation obstructions or a narrow width (horizontal features). Horizontal solid lines depict subreach boundaries and the dashed line represents ecoregion boundary. (b) Averaged daily fraction of flowing channel (FF ) in the extractable Hassayampa River (HR) reaches. Winter is defined as November to March and summer as July to September. (c) Average fraction of flowing days (FF ) in the HR along with flow thresholds in AZ/NM Mountains (ANM) and Sonoran Basin and Range (SBR) ecoregions. Annual flowing fractions less than the threshold indicate ephemeral status. In (b) and (c), raw data (gray lines) and smoothed results (black lines) are shown.

Spatial variations of streamflow presence averaged for the three WYs were compared to existing thresholds (Granato, 2010) for the two ecoregions in Figure 4c. The highly resolved determination of annual flowing days allowed a novel view of variations within and among subreaches. In the ANM, most of Reach 7B and 5 were ephemeral, with an average FF less than the threshold of 7.1% (Table S4 in Supporting Information S1). The fraction of reaches in the ANM ecoregion classified as intermittent was 9.8%. In contrast, subreaches in the SBR ecoregion had fewer flowing days, but a greater proportion of intermittent reaches due to the lower threshold of FF (3.6%). Large extents of Reaches 4, 3, and 1B exhibited a high number of flowing days and intermittent fractions (Table S4 in Supporting Information S1). Spatial variations in the characteristics of Reach 2A were identified with intermittent conditions observed upstream and downstream of a large riparian area known to be perennial (Katz et al., 2012). Further downstream, however, the flow regime became ephemeral, with the annual FF < 0.8% in Reaches 2B and 1A. Overall, the SBR ecoregion exhibited 16 km (or 12.9%) of intermittent rivers.

Discussion

Broader Applicability for IRES Management

Given that drylands constitute ∼40% of the Earth's land surface (Okin et al., 2015), the high spatiotemporal and global coverage from CubeSats offer the potential for determining streamflow presence where IRES have significant effects on hydrology, biogeochemistry, transport processes, and habitats. To demonstrate this broader applicability, we applied the method to 15 selected USGS gages in the SBR ecoregion outside of the HR (Figure S5 in Supporting Information S1). We found similar patterns of NIRDiff in most gages (12 out of 15, Figure S6 in Supporting Information S1), confirming the basic principles of the new method, except for sites where reaches were undetectable due to vegetation obstructions or a narrow width. With the optimal θ, the NIRDiff method performed well with respect to the fraction of flowing days (FF ), with an annual bias of 0.6% and a mean absolute error (MAE) of 11.0% when averaged across 12 sites (see details in Table S5 and Text S5 in Supporting Information S1). It is noteworthy that those gaging sites had larger values of θ than the gages in the HR, ranging from 0 to +0.098 (Table S5 in Supporting Information S1). A statistical analysis revealed a significant linear relationship (p < 0.0001) between F Gage and θ (Figure S7, Text S6 in Supporting Information S1), which provides means to estimate θ for a wide range of IRES. This relation explains why higher values of θ were found at SBR sites since FF is higher (62.4%) than at HR sites. This also suggests that θ can be obtained at a single gaging station and applied to reaches with approximately the same FF (Text S6 in Supporting Information S1) or extrapolated based on a linear regression established across a range of gaging sites with different FF (Figure S7 in Supporting Information S1).

Uncertainties and Limitations

As mentioned previously, it is important to carefully select the θ value before deriving the fraction of flowing days from PlanetScope imagery (F NIR). The positive linear relationship between optimal θ and F Gage suggests that using one regional θ for different rivers that have contrasting FF would lead to large errors in F NIR as compared to F Gage (see details in Text S6 in Supporting Information S1). Previous studies have used machine learning to evaluate natural or anthropogenic drivers of FF in non‐perennial rivers (e.g., Hammond et al., 2021; Messager et al., 2021). Such a method might be useful in future work to estimate the spatial distribution of θ across different river reaches using readily‐available ancillary data, for instance the aridity index (Messager et al., 2021). In addition, measurement uncertainties from PlanetScope imagery need to be carefully accounted for through quality control procedures prior to applying the streamflow detection method (see Text S7 in Supporting Information S1). Despite the near‐daily revisiting frequency, cloud contamination leads to missing data during the period when there might be flowing water. As a result, it is important to quantify the probability of flow presence when PlanetScope imagery is available (i.e., a subset of days in a year) and to evaluate if the record length is sufficient to characterize the streamflow regime.

Conclusions

By their nature, non‐perennial rivers are expected to have very high degrees of spatial and temporal variability in streamflow presence. However, direct evidence is often limited to: (a) time series at a few gaging stations and monitoring locations along a river (Gao et al., 2021; Gungle, 2005; Vivoni et al., 2006), and (b) mapping performed by field crews and citizen scientists during short campaigns (Allen et al., 2019; Lovill et al., 2018; Turner & Richter, 2011). In the Hassayampa River of central Arizona, for instance, the spatiotemporal variability of flow presence has previously been characterized by vegetation and other aquatic indicators on an annual basis at eight channel sites (Stromberg et al., 2009). While valuable, field methods are typically unable to detect rapid changes in water presence occurring over short distances due to infrequent storms. As shown here, CubeSat remote sensing can provide a comprehensive view of the spatiotemporal variability of streamflow presence over a large dryland river. In addition to estimating the annual flowing days with accuracy, derived products from the approach allow quantifying: (a) the fraction of a river that is flowing or has water presence for each day, and (b) the fraction of intermittent conditions within a river reach at 90 m resolution.

Robust streamflow detection with CubeSat remote sensing is feasible in dryland rivers due to rapid, high resolution detection of sharp spectral contrasts between water flow and surrounding dry areas. This novel approach can be readily used by scientists and environmental management agencies to determine streamflow regimes across global dryland regions. It complements field monitoring of IRES by helping to answer questions such as: “Where and how often is ground‐based sampling needed for streamflow determination?” This can improve the operational efficiency of field monitoring while reducing costs for streamflow assessments across large areas. Since the method is calibrated to available gaging data, and potentially with other field techniques, it builds confidence for its use in jurisdictional settings.

In recent years, there is a growing attention toward better understanding of streamflow permanence in the scientific and management communities (Jaeger et al., 2019, 2021; Mazor et al., 2021). The CubeSat detection approach presented here can be used as a novel source of data to evaluate and complement other efforts. For example, a recent review paper by Zimmer et al. (2020) pointed out that stream gages can produce zero‐flow readings despite the presence of water in the channel. CubeSat imagery can be used as an additional source of data to quantify the accuracy of zero‐flow readings and refine the spatial representation of streamflow permanence. Furthermore, the spatiotemporal products obtained from CubeSats can be used to calibrate and test physical and statistical models that simulate dryland streamflow regimes (e.g., Jaeger et al., 2014; Messager et al., 2021).

Supporting Information S1

Click here for additional data file.