Introduction

This report describes a continuation of work completed by the U.S. Geological Survey, in cooperation with the New York State Department of Environmental Conservation, and reported in Stagnitta and others (2024a). The low-streamflow statistics for the lowest annual 7-day and 30-day average streamflow that occurs (on average) once every 10 years (7Q10 and 30Q10) were calculated for 213 streamgages in unaltered watersheds and 79 streamgages in altered watersheds (hereafter referred to as “unaltered streamgages” and “altered streamgages,” respectively) across New York and hydrologically connected streams from bordering States. Stagnitta and others (2024a) defined altered streamgages as those with “documented alterations of observed daily streamflow due to human-related water use and management such as reservoir operations, surface water or groundwater withdrawals, diversions, engineered drainage systems, and impervious areas from urban development” (p. 2). The U.S. Geological Survey streamgage network in New York (https://waterdata.usgs.gov/ny/nwis/nwis) is well maintained and provides adequate spatial coverage for monitoring current and evaluating past streamflow conditions across the State (U.S. Geological Survey, 2024b). Historically, users of low-streamflow statistics, such as the New York State Department of Environmental Conservation (NYSDEC), have often required estimates of low-streamflow statistics in areas without available streamgages (Randall, 2010; Feaster and Lee, 2017; Randall and Freehafer, 2017; Lukasz, 2021). The NYSDEC uses low-streamflow statistics to determine pollutant loading limits for effluent from wastewater treatment facilities (NYSDEC, 1998). Currently, for facilities that are not near streamgages, low-streamflow statistics are estimated by methods such as the drainage area ratio between a correlated nearby streamgage with low-streamflow statistics available (Edward Schneider, New York State Department of Environmental Conservation, written commun., 2024). Low-streamflow statistics are also used for drought monitoring, watershed management, and dam operations (Smakhtin, 2001).

Development of multiple linear regression equations is a common approach to estimating low-streamflow statistics for ungaged basins (Feaster and Lee, 2017; Lukasz, 2021). Randall (2010) and Randall and Freehafer (2017) derived regression equations to estimate low-streamflow statistics at ungaged locations in New York; the data included in these studies were for the Susquehanna River Basin and the lower Hudson River Basin, respectively.

Furthermore, recent literature has shown that machine learning methods, such as random forest regression models, can provide improved model performance to estimate low-streamflow statistics (Eng and others, 2017; Worland and others, 2018; Eng and Wolock, 2022; DelSanto and others, 2023). The objective of this study was to develop random forest regression models that can be used to estimate the 7Q10 and 30Q10 flows for ungaged locations in New York (excluding Long Island) and to make the modeled output data available in StreamStats (U.S. Geological Survey, 2024a) to enable users to estimate low-streamflow statistics. Streamgages across Long Island were not used because its unique hydrogeology, urbanization, and regulation would require additional analysis beyond the scope of this study (Stagnitta and others, 2024a).

This report describes how the random forest regression models were developed, including the methods of basin delineation and removal of redundant or altered basins; the selection of basin characteristics as candidate predictor variables; and the process of training, tuning, and testing the models through machine learning methods. The report compares the modeled low-flow statistics with the observed statistics, presents the explanatory power of each basin characteristic in the 7Q10 and 30Q10 models, and shows the results of bias correction. Uncertainty, limitations, patterns, and interpretations of the results are also discussed. The low-flow statistics used to develop the models are available in Stagnitta and others (2024b) and discussed in Stagnitta and others (2024a). The basin characteristics developed for unaltered streamgage locations and R scripts (version 4.3.2; R Core Team, 2023) used to develop the models are available in Stagnitta and Woda (2025) and discussed in this report.

Study Area and Supporting Work

The scope of Stagnitta and others (2024a) included selection of streamgages in New York and hydrologically connected streams in bordering States, identification of unaltered and altered streamgages, data handling for missing and irregular annual minimum n-day values, a trend analysis to subset the annual minimum n-day time series for streamgages with significant trends, and calculation of the 7Q10 and 30Q10 statistics. The trend analysis performed in Stagnitta and others (2024a) identified streamgages with significant trends using a “Wilcoxon rank-sum hypothesis test to determine whether data from the most recent 30 years of record are statistically different from data from 30 years ago and earlier” (p. 5–6). Low-streamflow statistics were calculated for streamgages with a significant trend (p value less than or equal to 0.1) by using data from the most recent 30 years of record and for the remaining streamgages (p value more than 0.1) by using data from the entire period of record (Stagnitta and others, 2024a).

The study area for this study was the same as for Stagnitta and others (2024a), and the 7Q10 and 30Q10 statistics for the 213 unaltered streamgages were used to develop the models described in this report, for estimating low-streamflow statistics at ungaged locations (8 unaltered streamgages were identified as outliers and removed from model development, resulting in a total of 205 streamgages used for this study; see the “Outlier Streamflow Data” section). The 7Q10 statistics ranged from 0.03 to 945 cubic feet per second (ft³/s), and the median value was 5.21 ft³/s (Stagnitta and others, 2024b). The 30Q10 statistics ranged from 0.07 to 1,113 ft³/s, and the median value was 6.98 ft³/s (Stagnitta and others, 2024b). In addition, the drainage areas ranged from 0.6 to 7,782 square miles (mi²), and the median drainage area was 68 mi².

Methods

Random forest regression models were developed to estimate 7Q10 and 30Q10 statistics for ungaged locations using basin characteristics as predictor variables. A machine learning methodology was used to train, tune, test, and bias-correct the models, and a redundancy analysis was done to prevent data leaking from the training dataset to the testing dataset.

Redundancy Analysis

U.S. Geological Survey streamgages have a variety of purposes including, but not limited to, monitoring for flood or drought conditions, dam operations, and water supply management (Hester and others, 2006). Where monitoring interests or needs are greater, streamgages can be numerous within relatively small, localized areas of the State. One example is the approximately 1,600 mi² area upstream from the six New York City drinking water supply reservoirs in the Catskill Mountains region of New York, where 52 streamgages are presently active (June 2024).

Concentrated streamgaging networks can serve many important functions, but during development of a model in which two streamgages are nested or redundant (two nested streamgages represent the same stream reach, are close in proximity, and have similar drainage areas), there is a potential for data leakage, such that one nested streamgage is in the training dataset and the other nested streamgage is in the testing dataset (Meyer and others, 2018). Leakage of information from the training to the testing dataset produces an overfit model that does not accurately represent unseen data.

Veilleux and Wagner (2021) described a method to remove redundant basins to improve regional skew coefficients for a flood-frequency analysis. This method uses the standardized distance between centroids of two basins (SD) and the drainage area ratio of each basin (DAR) to compare two potentially redundant basins.

$SD=\frac{D_{ij}}{\sqrt{0.5\left(D{A}_i+D{A}_j\right)}}$

,

(1)

where

D_ij

is the distance between the centroid of the drainage basins upstream of the two streamgages i and j, in miles;

DA_i

is the drainage area of streamgage i, in square miles; and

DA_j

is the drainage area of streamgage j, in square miles.

$DAR=max\left[\frac{D{A}_i}{D{A}_j},\frac{D{A}_j}{D{A}_i}\right]$

,

(2)

where

DAR

is the maximum of the two values in brackets;

DA_i

is the drainage area of streamgage i, in square miles; and

DA_j

is the drainage are of streamgage j, in square miles.

Pairs of streamgages were flagged for potential redundancy if the SD was less than or equal to 0.5 and the DAR was less than or equal to 5. Figure 1 shows the delineated basins for the 213 unaltered streamgages and which basins were flagged as redundant. This analysis identified 69 streamgages with redundant basins, which were not included in the testing dataset or the cross-validation procedure to tune the model hyperparameters. The full dataset of 213 unaltered streamgages were considered when selecting the 70-percent training dataset to maximize the number of streamgages available to develop the models.

Model Development

Recent literature has shown that machine learning methods can outperform more traditional modeling approaches such as multiple linear regression for estimating low-streamflow statistics in ungaged basins (Eng and others, 2017; Worland and others, 2018; Eng and Wolock, 2022; DelSanto and others, 2023). For example, Eng and Wolock (2022) compared the performance of multiple linear regression with three machine learning methods (random forest, cubist regression, and support vector regression) to evaluate the performance of each method for several streamflow statistics across the contiguous United States. The random forest regression model estimated the 7Q10 and 30Q10 low-streamflow statistics more accurately compared to other methods. DelSanto and others (2023) focused their study on the Northeast region of the United States and compared multiple methods of estimating 7Q10, including multiple linear regression, random forest, neural networks, and a generalized additive model. Overall, the random forest regression method more accurately predicted measured 7Q10 values compared to other methods across a range of small (less than 15 mi²), medium (15–70 mi²), and large (more than 70 mi²) basins in the Northeast.

The model development for this study used R code available from Zipper and others (2021) and followed a workflow similar to the one they used to evaluate stream intermittency across several stream metrics using a random forest regression model. The R packages tidymodels (Kuhn and Wickham, 2020), ranger (Wright and Ziegler, 2017), and partykit (Hothorn and Zeileis, 2015) were used to train, tune, and test random forest regression models to estimate 7Q10 and 30Q10 statistics for unaltered streamgages in the study area.

Figure 2 shows the modeling workflow used to evaluate the unaltered streamgages (that is, to identify redundant basins and outlier streamgages; after their removal, a total of 205 streamgages were available for development of the model), the 70/30 training/testing dataset split, and the modeling procedure to train, tune, and test the model. The following steps explains a more detailed description of the modeling workflow.

1. 1. Log-transform the 7Q10 and 30Q10 dependent variables. Other transformations were explored such as dividing the low-streamflow statistics by the drainage area, but the log-transformation produced the most accurate output model results.

2. 2. Randomly split 70 (143 streamgages) and 30 (62 streamgages) percent of the dataset into training and testing datasets (streamgages with basins flagged as redundant are not included in the testing dataset but can be included in the training dataset).

3. 3. Train an initial random forest model by using a conditional random forest model to determine the most important predictor variables (described in more detail in the “Training—Variable Importance” section).

4. 4. Train a second random forest model, in which the number of predictor variables is incrementally increased, starting with the most important variables (from step 3), and determine the models with the lowest mean square error (MSE) and coefficient of determination (R^2; described in more detail in the “Training—Forward Feature Selection” section).

5. 5. Use the results from step 4 to tune the model hyperparameters—mtry, the maximum number of predictor variables in an individual tree; ntree, the maximum number of trees; and min_n, the minimum number of nodes or observations in an individual tree—using a tenfold cross-validation procedure (the resampled training dataset does not include redundant basins, which is described in more detail in the “Model Tuning” section).

6. 6. Test the tuned model with the 30-percent testing dataset and evaluate the performance of the model (described in more detail in the “Testing—Performance Metrics” section).

7. 7. Bias-correct the modeled output data by using an empirical cumulative distribution function to match modeled output with the observed data (described in more detail in the “Bias Correction” section).

Training—Variable Importance

Random forest regression models do not generally lose performance when linearly correlated predictor variables are included in the models (Konapala and Mishra, 2020). A random forest regression model is developed using hundreds of decision trees: a randomly selected subset of predictor variables is used within each decision tree, and the best predictors from the random subset are selected at each split (Ishwaran and others, 2010). Groups of correlated predictor variables that are important for estimating the dependent variable are split into separate trees to diminish the influence of multicollinearity but can ultimately be included as potential predictor variables in the final model (Ishwaran and others, 2010).

The initial models developed to calculate the importance of each predictor variable use a conditional variable importance parameter to select variables with a significance level (p value) of 0.05 or lower (the hyperparameter mincriterion equal to 0.95). The conditional random forest regression model measures the “mean decrease in accuracy” and accounts for both the main and interactive effects among predictor variables (Strobl and others, 2008). Also, the ntree hyperparameter for the cforest function was set to a high number (1,750 trees) to ensure that the same predictor variables were output with different randomly sampled training and testing datasets. After the initial models to determine the conditional variable importance for all 224 basin characteristics were developed for this study, 188 and 192 basin characteristics remained for the 7Q10 and 30Q10 models, respectively (basin characteristics less than or equal to 0.05 significance level remained).

Training—Forward Feature Selection

The forward feature selection method is a common approach in machine learning to reduce the total number of predictor variables used for a model (Khalid and others, 2014; Meyer and others, 2018). Removing predictor variables that may not be relevant or are potentially redundant improves model accuracy (Khalid and others, 2014). Several search methods exist for identifying the optimal models, such as forward and backward feature selection and compound, weighted, and random feature selection (Khalid and others, 2014). The forward feature selection method incrementally increases the number of predictor variables in successive model runs, and predictor variables for the model with the lowest error are used for additional model tuning and testing.

The methods for this study includes the following three steps:

1. 1. The output from the initial variable importance model was run to sort the predictor variables in order of importance.

2. 2. The individual models were run while the number of predictor variables within each model was increased incrementally.

3. 3. The model with the lowest MSE was determined.

The model with the lowest MSE is not always the best model to select, and models with fewer predictor variables are generally preferred (Helsel and others, 2020). Therefore, for this study, the total number of predictor variables within each model was considered. After forward feature selection models were developed and the total number of predictor variables was considered, the models with the 4th lowest (7Q10) and 5th lowest (30Q10) MSE, each of which had a total of 26 different predictor variables, were selected. A few models had lower MSE values (marginal improvements of log 0.01 ft³/s), but all these models had more than 30 predictor variables.

Model Tuning

The goal of model tuning is to provide an optimal tradeoff between underfitting and overfitting the dependent variable by tuning the model hyperparameters. A common approach to model tuning in machine learning applications is a k-fold cross-validation, in which the training dataset is randomly split into equally sized folds and each successive fold is held out as a validation dataset to tune the model (Meyer and others, 2018). In addition, Meyer and others (2018) determined that a “target-oriented” validation procedure, in which the spatial location of sites was considered, reduced overfitting and improved overall model accuracy. Sites that were in similar locations were separated from the random subset of the dataset for validating the models.

The random forest regression models included three hyperparameters to tune: the total number of predictor variables (mtry), the total number of trees (ntree), and the minimum number of observations or nodes within each tree (min_n). The models were tuned by a tenfold cross-validation procedure, in which 10 equally sized and randomly selected subsets of the 70-percent training dataset were held out individually to validate the model (redundant basins were removed from the 70-percent training dataset; a total of 80 streamgages remained). The cross-validation procedure allows for a grid search of the hyperparameters: mtry, which ranged from 1 to the maximum number of predictor variables determined from the forward feature selection for the 7Q10 and 30Q10 models (26 predictor variables for both 7Q10 and 30Q10), ntree, which ranged from 250 to 1,750, and min_n, which ranged from 3 to 25. The models with the lowest mean average error were selected, and the tuned hyperparameters were used for evaluating the performance of the models with the 30-percent testing dataset. The tuned 7Q10 model included a total of 8 predictor variables (mtry), 464 individual trees (ntree), and a minimum of 3 separate gaged observations within each tree (min_n). The tuned 30Q10 model included a total of 18 predictor variables (mtry), 464 individual trees (ntree), and a minimum of 3 separate gaged observations within each tree (min_n).

Testing—Performance Metrics

Performance metrics were used to evaluate how well each model fit to observed data from the testing dataset and included the following: the mean absolute error (MAE), the root mean square error (RMSE), the coefficient of determination (R²), and the Kling-Gupta efficiency (KGE). The KGE is similar to the Nash-Sutcliffe Efficiency in that the KGE accounts for bias, correlation, and variability as separate components (Gupta and others, 2009; Zipper and others, 2021). A KGE value more than or equal to −0.41 determines that a model is better than using the mean flow to estimate streamflows, and a perfect fit model has a KGE equal to 1 (Knoben and others, 2019 and Zipper and others, 2021).

There have been recent efforts to improve the interpretability of machine learning models, including the development of Shapley additive explanations (SHAP; Molnar, 2022). SHAP values are an estimation of Shapley values, a method derived from coalitional game theory, where features are “players,” and the prediction is the “payout” (Molnar, 2022). SHAP values are used to evaluate the explanatory power of the predictor variables on the individual predictions of the dependent variables (Lundberg and Lee, 2017). Several studies related to water resources have used SHAP values to aid in explaining the effects of predictor variables within a machine learning model related to water quality and quantity (Ransom and others, 2022; Madhushani and others, 2024). In this study, SHAP values were calculated and visualized by using the R packages fastshap (Greenwell, 2024) and shapvis (Mayer and Stando, 2024), respectively.

Results

The random forest regression models for finding the 7Q10 and 30Q10 statistics were evaluated by using a 30 percent testing dataset held out of the training and tuning procedures (Stagnitta and Woda, 2025). The streamflow values in the testing dataset ranged from 0.05 to 76.53 ft³/s for the 7Q10 model and from 0.08 to 91.99 ft³/s for the 30Q10 model. The ranges of the 7Q10, 30Q10, and basin characteristic values used for testing and training datasets for each model are shown in table 1.

Model Performance

The 7Q10 and 30Q10 models yielded evenly distributed homoscedastic residuals and an adequate fit to the observed data (figs. 3 and 4). The overall performance results for the 7Q10 model (MAE equal to 3.214 ft³/s, RMSE equal to 5.686 ft³/s, R² equal to 0.796, and KGE equal to 0.716) and the 30Q10 model (MAE equal to 3.677 ft³/s, RMSE equal to 6.133 ft³/s, R² equal to 0.853, and KGE equal to 0.811) confirmed that the models accurately predicted 7Q10 and 30Q10 values, and the prediction accuracies were similar to those from the random forest models used for estimating low-streamflow statistics, which are described in the literature (Eng and Wolock, 2022; DelSanto and others, 2023).

The models tend to overpredict at lower 7Q10 (less than about 2.7 ft³/s) and 30Q10 (less than about 1 ft³/s) observed values in log space (figs. 3A and 4A). In addition, the models tend to underpredict at higher 7Q10 (greater than about 50 ft³/s) and 30Q10 (greater than about 75 ft³/s) observed values (figs. 3B and 4B). Prediction intervals were determined by using quantile random forest regression to calculate the 5th and 95th percentiles of individual predictions from the model (Meinshausen, 2006; Wright and Ziegler, 2017); these intervals indicate the uncertainty in the model-predicted streamflows.

Figures 5 and 6 show, for 7Q10 and 30Q10, the mean absolute SHAP value and the range of individual SHAP values for basin characteristics (in log space), and the figures include both the training and testing datasets for the top 8 and 18 most important predictor variables for the 7Q10 and 30Q10 models, respectively. SHAP values represent the explanatory power of each feature (basin characteristic) on an individual 7Q10 or 30Q10 value prediction (Molnar, 2022). Larger average SHAP values indicate a larger influence of the predictor variable on the estimate of the low-streamflow statistics.

The basin geometry characteristics Drainage_Area, Basin_Stream_Length, Basin_Perimeter, and Longest_Flow_Path were the most influential predictor variables for both the 7Q10 and 30Q10 models. These variables were well distributed from low to high SHAP values, indicating these variables are important for predicting both large and small 7Q10 and 30Q10 values.

Several additional basin characteristics differed in prediction influence for 7Q10 and 30Q10 models. Three basin characteristics were significant only for predicting 7Q10 values:

• • Mean_Oct_Tmax—basinwide average of maximum October temperatures,

• • Slope_LFP_Ratio—the ratio of the basin slope to the longest flow path, and

• • Mean_Jun_Tmax—basinwide average of maximum June temperatures.

Thirteen basin characteristics were significant only for predicting 30Q10 values:

• • SSURGO_Hydric_Class_Presence—the basinwide percentage of soils from the hydric class presence category (an indicator of wetland soils),

• • STD_Apr_Tmin—the basinwide standard deviation of minimum April temperatures,

• • Mean_Jun_Tmean—basinwide average of mean June temperatures,

• • PS_Coarse_Loamy—basinwide percentage of soils classified as coarse-loamy,

• • Mean_Mar_Precip—basinwide average of March precipitation,

• • Mean_Apr_Tmean—basinwide average of mean April temperatures,

• • STD_Oct_Tmax—the basinwide standard deviation of maximum October temperatures,

• • STD_Aug_Tmax—basinwide standard deviation of maximum August temperatures,

• • STD_Apr_Tmax—basinwide standard deviation of maximum April temperatures,

• • Surficial_LF_Alluvial—basinwide percentage of the basin classified as alluvial surficial geologic materials,

• • STD_Jun_Tmin—basinwide standard deviation of minimum June temperatures,

• • STD_May_Tmin—basinwide standard deviation of minimum May temperatures, and

• • Slope_LFP_Upper—the slope between the 0- and 50-percentile elevations along the longest flow path.

One basin characteristic was significant for predicting both 7Q10 and 30Q10 values:

• • SSURGO_HSGA—the basinwide percentage of soils from hydrologic soil group A (well-draining soils).

Bias Correction

Machine learning regression models, like many other related models, can systematically overestimate small values and underestimate large values (Belitz and Stackelberg, 2021). Belitz and Stackelberg (2021) evaluated six methods for correcting bias in predictions from machine learning regression models and determined that a simple empirical matching of the distributions of predicted values to observed values most accurately predicted hold-out datasets of all methods evaluated. Here, the training dataset was used to match the observed empirical distribution to the predicted values, and this distribution was then used to bias-correct predicted values from the testing dataset. Figures 7 and 8 show the cumulative distribution for the observed, modeled, and bias-corrected training and testing datasets. Figures 7A and 8A show that the modeled training dataset fits the observed distribution well before bias correction, whereas figures 7B and 8B show that when the models are used for testing with unseen data, the bias in predictions is pronounced. The testing dataset shows the bias in machine learning regression models to overestimate small values and underestimate large values, and the bias correction improves the predictions of these extreme values (figs. 7B and 8B). In addition, a bias correction in log space ensures that large 7Q10 or 30Q10 values do not dominate the distributions (Belitz and Stackelberg, 2021).

Discussion

Figures 3 and 4 show the 95-percent prediction interval for each modeled prediction of 7Q10 and 30Q10, and the potential ranges of predicted values were large (Stagnitta and Woda, 2025). The bias correction described in the previous section can help to reduce the error in predictions for new, unseen datasets, but there is still a level of uncertainty in the modeled output data. The percentages of observed values from the testing dataset outside the range of the prediction intervals were 4.8 (3 out of 62 streamgages) and 8.1 (5 out of 62 streamgages) for the 7Q10 and 30Q10 models, respectively.

The models were developed using a total of 205 unaltered or minimally altered streamgages. The models are not intended to be used in locations where there are potential human-induced effects on streamflow, such as withdrawals, dam operations, canal diversions, or high urbanization. If a user predicts a low-streamflow statistic at an ungaged location where there is a known flow alteration, the user should account for the potential effects of alteration.

Random forest models can accurately model low-streamflow statistics with smaller datasets (a few hundred observations), but model accuracy improves as the amount of observed data available to train the model increases. Biau (2012), found that model accuracy generally improved when at least 500 observations were included while training the model. As additional streamgages come online and newer data become available for model training, further updates and improvements to the models developed for this study may be warranted.

The models developed for this study were created to prioritize the prediction of 7Q10 and 30Q10 for ungaged locations across New York. These models were not developed to be used in basins with drainage areas outside the ranges in the training dataset (table 1). In addition, the models are not intended to be used to identify the physical drivers influencing low streamflow (correlated predictor variables in the final models diminish the usability of the models to describe the drivers). Because some variables in the final models were correlated, further study of the influence of predictor variables (described below) could improve understanding of the drivers of low-streamflow statistics in the study area.

The basin characteristics Drainage_Area, Basin_Stream_Length, Basin_Perimeter, and Longest_Flow_Path were the four most important predictor variables for the 7Q10 and 30Q10 models. These basin characteristics describe the physical topography and the overall length of streams within the basin and are the main physical drivers that determine the amount of runoff from precipitation that could be feasibly received within a network of streams. The models showed that larger values for these four predictor variables had larger low-streamflow statistics and that smaller variables had lower low-streamflow statistics.

In addition, the 7Q10 model included the slope between the most upstream and downstream points of the longest flow path divided by the mean basin slope (Slope_LFP_Ratio; determined by the intersection of basin polygon coverages and a digital elevation model). The 30Q10 model included the slope between the 0- and 50-percentile elevations of the longest flow path (Slope_LFP_Upper). The 7Q10 model showed that basins with a higher Slope_LFP_Ratio (basins in which the slope of the longest flow path and the mean basin slope were similar) had higher low-streamflow statistics. Generally, these were larger basins, in which the mean basin slope is less influenced by extreme elevations within the basin than basins with higher Slope_LFP_Ratio values (these were generally smaller basins). The Slope_LFP_Upper showed a similar pattern to the Slope_LFP_Ratio but to a much lesser extent: lower values of Slope_LFP_Upper corresponded to higher low-streamflow statistics.

Several basin characteristics related to temperature were included for the 7Q10 (Mean_Oct_Tmax and STD_Apr_Tmin) and 30Q10 (STD_Apr_Tmin, Mean_Jun_Tmean, Mean_Apr_Tmean, STD_Oct_Tmax, STD_Apr_Tmax, STD_Aug_Tmax, STD_May_Tmin, and STD_June_Tmin) models, respectively. The 7Q10 model showed that basins with higher maximum October temperatures tended to have lower 7Q10 values, indicating the influence of high temperatures in reducing streamflow during earlier fall months (September and October). The 30Q10 model showed a similar influence, to a lesser extent, of spring and early summer temperatures on streamflow. The 7Q10 and 30Q10 models showed that temperatures with a higher standard deviation or variability across the basin (on the basis of PRISM grid cells) of temperatures in spring, summer, and fall months for 7Q10 and 30Q10 models tended to have higher low-streamflow values. The higher variability temperature values occurred in basins that were largely forested and included both low and high elevations within the basins (Catskill and Adirondack Mountain regions). This may indicate basins with potentially more groundwater storage could provide more sustained low-streamflow values throughout the year. The 30Q10 model included the Mean_Mar_Precip, where lower values indicated lower low-streamflow statistics, which may indicate the importance of late winter or early spring precipitation in recharging groundwater-dominated basins.

In addition, basin characteristics related to soils and surficial geology were included with the 7Q10 model (SSURGO_HSGA) and 30Q10 model (SSURGO_Hydric_Class_Presence, SSURGO_HSGA, PS_Coarse_Loamy, Surficial_LF_Alluvial). The SSURGO_HSGA variable shows the percentage of the basin covered by hydrologic soil group A and is defined by well-draining soils. Basins with a low percentage of well-draining soils tended to have relatively low low-streamflow statistics. High SSURGO_HSGA percentages could be indicative of basins less able than others to infiltrate runoff to groundwater, and in such basins, groundwater would have less influence on streams during the seasonal low-streamflow period.

The 30Q10 model included more basin characteristics related to soils and surficial geology than the 7Q10 model. The 30Q10 low-streamflow statistics were calculated from a 30-day moving average, which is a longer period than the 7-day moving average used for the 7Q10 statistics. The soils and surficial geologic basin characteristics indicate a basin’s ability to store and mediate the flow of groundwater, which is likely more applicable to low-streamflow statistics calculated from longer moving average periods than from shorter periods. The percentage of the basin covered by the surficial coarse and loamy geologic soil class (PS_Coarse_Loamy) showed a pattern similar to that of the SSURGO_HSGA percentage but to a lesser degree: lower values corresponded to lower low-streamflow statistics. The percentage of the basin classified as “hydric” soils (SSURGO_Hydric_Class_Presence) is an indicator of wet soils and could also indicate the presence of wetlands; basins with higher percentages of hydric soils tended to have lower low-streamflow statistics than other basins. Generally, basins with wet soils and high percentages of wetlands store more water and could limit the amount of water available during the low-streamflow period. The percentage of the basin covered by the surficial geologic class of alluvial materials (Surficial_LF_Alluvial) showed a pattern similar to that of the SSURGO_Hydric_Class_Presence percentage but to a lesser degree: higher values of alluvium coverage corresponded to lower low-streamflow statistics.

The final models developed generally were well fit to the observed data, with an R² of 0.796 for the 7Q10 model and an R² of 0.853 for the 30Q10 model. The random forest regression model results are consistent with other studies across the Northeast and conterminous United States, which show that random forest regression models have higher prediction accuracies than models developed by linear regression methods (Eng and Wolock, 2022; DelSanto and others, 2023). DelSanto and others (2023) developed a random forest regression model for predicting 7Q10 across 106 streamgage locations from the Northeast, and they reported an R² of 0.61 from a leave-one-out cross-validation procedure.

Additionally in this study, a bias-correction method was developed to remove bias by matching the distribution of predicted to observed values. The matched distribution can then be used to bias-correct new model results for ungaged locations. The bias correction will automatically be applied to output model results available in StreamStats.

StreamStats Web Application for Modeled Results in Ungaged Locations

The U.S. Geological Survey developed and maintains the interactive web application StreamStats (U.S. Geological Survey, 2024a). StreamStats provides functionality for users to delineate basins and obtain peak- and low-streamflow statistics for gaged and ungaged locations across the United States. The at-site 7Q10 and 30Q10 low-streamflow statistics calculated by Stagnitta and others (2024a) were the first low-streamflow statistics that were made available in StreamStats for New York. The interactive web application includes functionality to obtain predictions of 7Q10 and 30Q10 by using the models developed for this study in ungaged locations across the State. The tool includes functionality that allows users to click on a stream grid cell to delineate a basin and to then obtain basin characteristics and predictions for 7Q10 and 30Q10 statistics for ungaged locations. For more information about StreamStats and how to use the tool, refer to the help section of the web interface (U.S. Geological Survey, 2024a).

Summary

Low-streamflow statistics are often needed for applications related to permitting and drought monitoring in locations without available streamgages. The U.S. Geological Survey, in cooperation with the New York State Department of Environmental Conservation, used basin characteristics as predictor variables to develop random forest regression models to predict low-streamflow statistics at ungaged locations in New York, excluding Long Island but including hydrologically connected basins from bordering States. The models were trained, tuned, and tested to ensure that they performed well across the State. A redundancy analysis was performed for all basins with unaltered streamgages to flag those that were close in proximity and had similar drainage areas. The redundancy analysis flagged 69 gaged basins, which were removed from the testing dataset and cross-validation procedures to prevent data leaking from the training dataset to the testing dataset. In addition, 8 streamgages were flagged as outliers, and a total of 205 streamgages were available for model development.

Basin characteristics were developed across the following categories: basin geometry, climate, land cover, soils and surficial geology, and other characteristics. The basin characteristics were used as predictor variables in the model. The most important predictor variables were determined by developing a model to calculate initial conditional variable importance, and the total number of predictor variables was further reduced by using a forward feature selection method. The physical basin characteristics related to the drainage area, total stream length, the perimeter of the basin, and the longest flow path were the top predictor variables for the models of the lowest annual 7-day and 30-day average streamflow that occurs (on average) once every 10 years (7Q10 and 30Q10). Additional basin characteristics related to slope, temperature, precipitation, soils, and surficial geology were also important predictors for the models.

The models were trained to determine importance values of the predictor variables and to reduce the number of predictor variables for tuning by using a forward feature selection method. The models were tuned by using a tenfold cross-validation procedure to determine the optimal maximum number of predictor variables within each tree (mtry), the total number of decision trees within the random forest model (ntree), and the minimum number of observed data points within each tree (min_n). The models were evaluated by using the testing dataset held out from training, and results indicted well-fit generalized models, with coefficients of determination (R²) of 0.796 and 0.853 for the 7Q10 and 30Q10 models, respectively. The output model results for ungaged locations within the study area are available in StreamStats.