We calculated observed annual values of 24 functional flow metrics (FFMs) describing 5 functional flow components (fall pulse flows, wet season baseflows, wet season peak flows, spring recession flows, and dry season baseflows) from reference gage records in California (Figure 1). We then characterized the watershed above each reference gage using a suite of physical and climatic variables from publicly available data sources. Next, we used a machine learning approach to relate the watershed variables to functional flow metrics, developing a total of 24 models (one for each functional flow metric). The predictive performance of each model was then evaluated by comparing predictions of functional flow metrics with observations at gages excluded from model training. Finally, we used the models to predict the natural range of values of each FFM at all stream reaches in California’s stream network, using the same set of predictor variables calculated for the catchment of each stream reach. The details of each step are provided below.

All gages operated by the U.S. Geological Survey (USGS) in California were screened to identify those considered to be reference-quality, following methods described by Zimmerman et al. (2018). Briefly, the watershed above each gage was evaluated using GIS-based methods and visual inspection of aerial imagery to exclude sites with evidence of significant human activities, including water diversions and storage reservoirs, intensive agriculture and forestry practices, dense road networks, and extensive impervious surfaces. We also reviewed USGS published annual data reports for each gage that note the influence of significant anthropogenic activities on observed flow records (Falcone et al., 2010). Through a subsequent manual screening process, several gages were removed from the analysis that exhibited irregular, impaired, and or aseasonal flow patterns. In total, we identified 216 reference gages in California, which included both active stations located on relatively pristine streams and gages with historical observations that pre-dated significant anthropogenic disturbances (e.g., prior to dam construction). We included an additional 3 gages located below dams for which reconstructed unimpaired flow data were available in order to increase the physiographic range represented in the dataset (California Department of Water Resources [DWR], 2007), bringing the total number of reference gages to 219 (Figure 1; Supplementary Table S1).

The resulting reference gage set includes periods of record as early as 1950 and as recent as 2015, with an average period of record of 33 years and ranging from 6 to 65 years (Supplementary Table S1). These gages are well distributed across the diversity of river types in California, including snowmelt (n = 25), rain (n = 125), and mixed snow-and-rain (n = 69) hydrologic regimes, following a simplified version of a stream classification scheme developed by Lane et al. (2017). The gages are located on streams with drainage areas ranging from 5 to 9,340 square kilometers and are distributed throughout California, with the exception of the arid southeastern corner of the state (Figure 1), where most streams are ephemeral and no reference gages are present. We confirmed that reference gages located in close proximity were separated by intervening tributaries or had distinct periods of record.

Complete years of daily flow data for all reference gages were downloaded from the USGS National Water Information System (USGS, 2017). Annual FFMs were then calculated from daily flow records using signal processing algorithms designed to characterize seasonal flow features of the annual hydrograph. The approach to calculate annual timing metrics detailed by Patterson et al. (2020) is as follows: A high standard deviation Gaussian filter was applied to daily streamflow time series to detect dominant peaks and valleys from the annual hydrograph. Localized search windows were set around hydrologic features of interest (e.g., annual peak flow). A low standard deviation Gaussian filter was then applied to the observed daily flow in the search window to identify seasonal shifts in the hydrograph, based on slope breaks in the derivative of a fitted spline curve. Break points were used to quantify the timing metrics for the wet season, dry season, and spring recession periods, from which seasonal magnitude, duration, frequency, and rate of change metrics could then be calculated. The spring recession rate was calculated as the median daily rate of change in flow from the start date of the spring recession until the start of the dry season, considering only days with negative change to omit storm events during the recession period. Peak flow magnitudes were calculated as the long-term annual flood exceedance flow associated with the 2-, 5-, and 10-years recurrence intervals. We also calculated peak flow duration (cumulative number of days in which this peak flow magnitude is exceeded) and frequency (number of times the flow magnitude is exceeded), in each of the years in which a flood of a given recurrence interval occurred. Following these methods, we calculated the values of 24 FFMs, describing 5 functional flow components, at each reference gage (Table 1; Figure 1).

Next, we used a GIS-based approach to calculate over 150 variables related to physical attributes of the watershed above each reference gage using publicly available geospatial datasets (Supplementary Table S1). These included variables related to topography (e.g., elevation, slope, and aspect, etc.), dominant geology and soil types (e.g., granitic, volcanic, or sedimentary, and mean content of clay, sand, and silt, etc.), and watershed hydraulic properties (e.g., topographic wetness index, baseflow index, mean depth to water table, etc.). We also included time-varying climatic variables including mean monthly temperature and precipitation from the 800-m PRISM dataset from Daly et al. (2008), as well as expected monthly runoff from McCabe and Wolock (2011). These climate variables were expressed as monthly, seasonal, and annual values for each year of FFM observations at a reference site, as well as for years preceding the FFM observations (Supplementary Table S2).

Random forest (RF) models (Cutler et al., 2012) were developed for each FFM. For most FFMs, observed values were calculated for each year of the reference period of each gage. For peak flow magnitude FFMs, single values for the 2-, 5-, and 10-years recurrence interval flood were estimated at each gage. Each RF model specified a FFM as the response variable and a total of 182 watershed and climate variables as predictor variables (Supplementary Table S2; Carlisle, 2022). All models were run using 2000 trees and default parameters with the randomForest function in the randomForest package, version 4.6 (Liaw and Wiener, 2018) in R (R Core Team, 2020).

Random forest models include a resampling routine that provides estimates of model performance comparable to what is obtained from independent validation data. However, because our dataset included repeated observations of FFM values from each of the 219 reference sites, the replicate datasets generated from RF’s internal sampling could produce overly optimistic estimates of model performance. We therefore used a leave-one-out cross-validation approach to estimate model performance, in which each reference site (including observations for all years of record) was excluded in turn from a calibration dataset, following methods by Eng et al. (2017); Zimmerman et al. (2018). The trained model was subsequently used to predict FFM values at the excluded reference site. We retained the 10th, 25th, 50th, 75th, and 90th percentiles of the predictions generated by the 2000 trees for each excluded reference site. For each model iteration, we also identified the most influential predictor variables, based on their Gini index (Cutler et al., 2007), which measures the loss of model predictive accuracy when that variable is excluded. Higher values indicate greater importance in contributing to the accuracy of the models.

To assess model performance, we compared predicted FFM values with observations at sites excluded from model training. We restricted the assessment to sites with 20 or more observations (i.e., 20 years of record) and calculated several model performance criteria to limit the risk of flawed interpretation resulting from the use of a single performance metric (Clark et al., 2021). We calculated performance criteria that provided measures of both the dispersion and central tendency of model predictions in comparison to observed values. First, we compared the distribution of observed to predicted values of FFMs by calculating the percent of annual observed values at a site that fell within the predicted interquartile range (IQR, range between the 25th to 75th percentile values) and the inter-80th percentile range (I80R, range between the 10th to 90th percentile values) for that site. The mean of these percentage values across all sites was used to assess the overall degree to which the distribution of observations aligned with the predicted range of each metric. Models with perfect performance would have percentage values of 50% for the IQR criterion and 80% for the I80R criterion, indicating that, on average, 50% and 80% of the observed values fall within the predicted IQR and I80R, respectively. Models that under-estimate the natural range of variation in FFMs would have values below 50% and 80%, respectively, and models that over-estimate the range of variation would exceed these values.

To evaluate accuracy for the central tendency of model predictions, we also compared the median value of observations to the median value of predictions at each site. The paired values were used to calculate several “goodness-of-fit” criteria commonly used in hydrologic model performance assessment (Moriasi et al., 2007; Eng et al., 2017): the observed-to-expected ratio (O/E), the coefficient of determination (r²), percent bias, and Nash-Sutcliffe Efficiency (NSE). We then calculated the mean value of each performance criterion across all sites. For the peak flow magnitude metrics, we only calculated the performance measures of central tendency because only single values were available for each site (i.e., 2-, 5-, and 10-years recurrence interval peak flows). Due to the skewed distribution of peak flow frequency and duration metrics to low values, measures of central tendency were unreliable. For those metrics, we excluded observations with zero values and considered only the distribution of observations relative to the predicted range of values, by calculating the percentage of observations falling within the predicted IQR and I80R, as described above.

To evaluate model performance across all criteria, we standardized the values of all calculated criteria between 0 (poor performance) and 1 (perfect performance). To scale O/E values, we retained values less than 1 and calculated the inverse of those greater than 1. To scale percent bias, we subtracted values from 100 and then divided by 100. NSE values less than 0 were set to 0 and no changes were made to the r² values. To scale the IQR criterion, the absolute value of difference between the calculated value and 50 was divided by 50 and subtracted from 1. Similarly for I80R, the absolute value of difference between the calculated value and 80 was divided by 80 and subtracted from 1. We then developed a composite performance index by averaging the values of all six criteria. We assigned a qualitative performance rating to the composite performance index values excellent (>0.9), very good (0.81–0.9), good (0.65–0.8), satisfactory (0.5–0.64), and poor (<0.5) model performance, following guidelines similar to Moriasi et al. (2007).

Finally, we evaluated spatial bias in model performance by separating reference gages into stream classes (Lane et al., 2017). We grouped gages into one of three classes based on their dominant hydrologic characteristics: snowmelt, rain, and mixed snow-and-rain. We then compared model predictions and observed data from reference gages occurring within each stream class, using the same set of performance criteria, and again calculated the composite performance index for each metric.

After the model performance evaluation, we used the RF models to predict the natural range of functional flows for all stream reaches in California. We trained final models (n = 2000 trees) with the full set of reference gages to include the maximum amount of information possible. We then calculated the same set of watershed and climate variables used in model training, obtained from Wieczorek et al. (2018), for 142,509 natural stream reaches (mean length = 2.1 km; sd length = 2.0 km) represented by the National Hydrography Dataset for California (NHDPlus, Version2) (Horizon Systems Corporation, 2012). These data were used to predict FFM values from 1950 to 2015 at each stream reach from the trained RF models. The 10th, 50th, and 90th percentiles of model predictions for each FFM were calculated for each stream reach. Predicted ranges were compiled for all years (1950–2015) and for all dry, moderate, and wet water years and made available on a public website (California Environmental Flows Working Group [CWFWG], 2021). Reported values represent the expected natural range of FFMs at each stream, also accounting for model prediction uncertainty.

Climate variables were generally the most influential predictors in the FFM models, although physical catchment variables were important for some metrics (Table 2; Supplementary Table S3). For the fall pulse metrics, climate variables including precipitation, temperature, and runoff for fall season months (e.g., Oct, Nov) were consistently among the most influential variables. The rainfall-runoff erosivity index, which reflects the estimated amount and rate of runoff produced by a storm (Renard et al., 1997), was also influential in predicting fall pulse flow duration. Wet season baseflow metrics were strongly influenced by monthly climate variables corresponding to winter months (e.g., December runoff) and multi-annual antecedent precipitation, runoff, and temperature variables were generally the most important variables in the models. Precipitation and runoff had a stronger influence on wet season baseflow magnitudes and duration, whereas temperature had a stronger influence on wet season baseflow timing.

Peak flow magnitudes, including 2-, 5-, and 10-years recurrence interval peak flow metrics, were most influenced by the catchment’s long-term mean annual runoff (Gebert et al., 1987) as well as mean maximum and mean annual precipitation (Table 2; Supplementary Table S3). Peak flow duration–the number of days in a year in which flows exceeded a peak flow threshold–was most influenced by monthly precipitation variables in the winter months and the previous water year, catchment mean elevation, and the hydrologic landscape region in which the catchment predominately occurs (Wolock, 2003a). Peak flow frequency—the number of peak flow events in a year of a given recurrence interval—was also most influenced by precipitation in winter months and antecedent year, as well as the catchment’s rainfall-runoff erosivity index and groundwater recharge index (Wolock, 2003b).

Annual precipitation and runoff had the greatest influence on spring recession magnitude, whereas mean temperatures in the spring months and spring season had the greatest influence on spring recession timing (Table 2; Supplementary Table S3). The duration of the spring recession flow period was most influenced by catchment elevation, the clay content of catchment soils, winter precipitation, and monthly runoff in December and June, near the start, and end of the wet season, respectively and the rate-of-change by mean runoff observed over the most recent four-year period. Spring rate-of-change was most influenced by runoff variables, including annual runoff for the water year and multi-annual antecedent periods.

Runoff variables were also influential in predicting dry season baseflow characteristics. Runoff in the summer months and seasonally-averaged runoff were both important in predicting dry season baseflow magnitudes. The catchment groundwater recharge index was also influential for the dry season high baseflow metric. Catchment variables were most important in predicting the timing of the dry season. Influential variables included the mean and max catchment elevation, the baseflow index, and soil properties (Supplementary Table S3). Similar to the fall and spring duration metrics, duration of the dry season was most influenced by catchment properties, including elevation, spring precipitation, and the catchment erodibility index (K-factor).

Overall, the FFM models had high predictive accuracy, with all 24 metrics exhibiting excellent (composite performance index [CPI] > 0.9 for 7 metrics), very good (0.8 < CPI ≤0.9 for 12 metrics) or good (0.65 < CPI ≤0.8 for 5 metrics) performance (Figure 2; Supplementary Table S4). The models performed well in predicting fall pulse flows, including the magnitude (CPI = 0.85), timing (CPI = 0.80), and duration (CPI = 0.70). The slightly lower CPI for fall pulse duration was driven by low NSE and r² performance criteria values (<0.25). This was the result of the limited range of whole number values in the observation record (median fall pulse duration of 2–7 days among all sites), such that slight deviation of predictions (i.e., 1 or 2 days) caused NSE and r² values to substantially decrease.

The models for wet season baseflow accurately predicted all metrics, especially the wet season low magnitude metric (CPI = 0.91), and exhibited excellent performance (CPI >0.9) in predicting peak flow magnitudes for 2-, 5-, and 10-years flood recurrence intervals. Model performance for within-year flood frequency and duration were also considered good, very good, or excellent. The model for the 10-years flood frequency tended to overestimate the observed range of variation (i.e., a higher proportion of observed values fell within the predicted interquartile range than expected; Supplementary Table S4), although overall model performance was still good. Model performance was very good (CPI >0.8) for spring recession flow magnitude, timing, duration, and rate-of-change. The models were very good or excellent in predicting all dry season baseflow metrics, including the median (CPI = 0.92) and high baseflow magnitudes (CPI = 0.92), dry season timing (CPI = 0.90), and dry season duration (CPI = 0.83).

When model performance was assessed by stream class, the CPI deviated from those obtained when all streams were evaluated together (Figure 2). For snowmelt-dominated streams, the models performed less well in predicting timing and duration metrics, including for the fall pulse, wet season, spring recession, and dry season. However, only the fall pulse timing model was considered “poor” performing for the snowmelt stream class. Model performance declined for some metrics in the mixed snow-and-rain and rainfall-dominated classes, but all models were considered at least satisfactory and most were very good or excellent (Figure 2). Stream gage records were insufficient to evaluate the performance of the 10-years flood duration and frequency metrics by stream class.

Based on the overall satisfactory performance of the FFM models, predictions of expected, natural FFM values were generated for all stream reaches in California using models calibrated with the full set of reference gages. Model predictions were compiled in a geospatial database and made available through an online mapping tool to allow users to visualize and download estimates of natural FFM values for any stream reach in the state CEFWG, 2021 (Figure 3).