Modeling Water Yield : Assessing the Role of Site and Region-Specific Attributes in Determining Model Performance of the InVEST Seasonal Water Yield Model

Simple hydrological models, such as the Seasonal Water Yield Model developed by the Natural Capital Project (InVEST SWYM), are attractive as data requirements are relatively easy to satisfy. However, simple models may produce unrealistic results when the underlying hydrological processes are inadequately described. We used the variation in performance of the InVEST SWYM across watersheds to identify correlates of poorly modeled outcomes of InVEST SWYM. We grouped 749 watersheds from across North America into five bioclimatic regions using nine environmental variables. For each region, we compared the predicted flow patterns to actual flow conditions over a 15-year period. The correlation between the modeled and actual flows was highly dispersed and relatively poor, with 92% of r2 values less than 0.5 and 42% less than 0.1. We linked cryospheric variables to model performance in the bioclimatic region with the poorest model performance (the Low elevation Boreal Sub-humid region—LeBSh). After incorporating cryospheric conditions into the InVEST SWYM, predictions improved significantly in 30% of the LeBSh watersheds. We provide a relatively straightforward approach for identifying processes that simple hydrological models may not consider or which need further attention or refinement.


Introduction
Globally, human populations are dependent upon water resources for their life-sustaining supply of fresh water as well as for their economic value [1].The ability to reliably model watershed processes is a critical component in allowing decision-makers and managers to maintain and improve aquatic ecosystem services (e.g., drinking water supply, recreational opportunities, and energy production) [2].The use of such models is also beneficial for predicting potential outcomes resulting from changes in precipitation or extreme climate events.For example, models can be used to identify areas susceptible to flooding or drought [3] under changing climate scenarios.However, hydrological processes are difficult to model as they are affected by complex factors such as climatic conditions, soil characteristics, land cover type, and topographic features, which vary spatially and temporally [4,5].In the last 20 years, the geographic information system technology has allowed the development of physical hydrological models to enable the inclusion of these complex factors [6] such as Penn State Integrated Hydrologic Modeling System Soil and Water (PHIM), the Soil and Water Assessment Tool (SWAT) [7] and the Soil and Water Integrated Model (SWIM) [8].Models can perform well when the appropriate data is available.Nevertheless, in the absence of robust and readily available datasets, simpler models are attractive to guide management and mitigation efforts aimed at optimizing ecosystem services.
All hydrological models present uncertainty related to their modeled processes and underlying structure [9][10][11].While the performance of process-based models should be robust against changing conditions, model structures can become invalid if the dominant processes underlying the model change [10,11].As above, the data requirements needed to run some of these models, and obtain useful results, are sometimes complex and/or difficult to satisfy.These requirements can make them difficult or even impossible to use for some regions, primarily due to the lack of available data for the region being studied.Therefore, before applying a hydrological model, it is necessary to verify the transferability of the model structures under different environmental conditions [11][12][13] and to determine which inputs have the greatest effect on the model results for each region.This verification would allow users to focus their efforts on collecting key data of sufficient quality to obtain meaningful results.
In this context, simpler models are potentially more attractive due to their ease of use and reduced data requirements [14].For example, a suite of tools developed by the Natural Capital Project as the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) framework [15] have become a popular choice in the analysis of a diverse set of ecosystem services.Ochoa and Cardona [16] found that, since 2009, InVEST has been the most used tool for modeling ecosystem services.The InVEST Seasonal Water Yield Model (SWYM) component has been used to asses water yield [16] in geographically diverse locations such as Hawaii, Indonesia, and China [6,17,18].The use of InVEST SWYM allows the exploration of different scenarios such as climate and/or land use change [19,20] as do other complex models such as SWAT.However, complex models need post-processing to allow decision-makers to assess the ecosystem services studied.In contrast, the InVEST SWYM model produces results that ultimately refer to the ecosystem services themselves and present them graphically for ease of visualization, permitting decision-makers to focus on the outcomes quickly and easily [6,14].Nonetheless, the minimal data requirements of a simpler model, particularly in combination with poor quality data, may produce results that do not represent reality, demanding extensive and complicated calibration and verification processes that may or may not improve results significantly.This result would ultimately defeat the purpose of using an easily applied model and is particularly important when applying the same model structure across regions with different environmental conditions.
For the purposes of our study, we focused on the InVEST SWYM as this model allows users to rank specific parcels (or pixels) of land for conservation or development based on their relative contributions to specific components of the hydrological cycle [15].That is, the InVEST SWYM model partitions precipitation into either quick flow or base flow (runoff versus groundwater recharge) by calculating a water balance for each individual pixel of the watershed of interest [14,21].The information required by the SWYM is easily obtained globally from publicly available data sources and includes monthly precipitation, topography, evapotranspiration, land-use, soil type, and land-cover data.Our objective was to analyze the performance of the InVEST SWYM across several environmental regions and use differences in performance across regions to isolate key variables that would improve the model.We propose that the results will differ depending on regional differences in environmental features, which will be especially critical in widely applied models such as InVEST SWYM due to the simplification of their processes, variables, and the parameters they incorporate in order to facilitate implementation.To test our hypothesis, we grouped 749 watersheds from across North America into five groups according to nine environmental variables and analyzed (1) how well the river flows modeled with InVEST SWYM corresponded with historical data in a regional context, and (2) assessed which environmental variables had the greatest influence on the modeled results.

Materials and Methods
Our methodological approach followed three main steps, which we detail below and outline in Figure 1.
Figure 1.Flowchart of methods.The main steps of the methodology (dark grey boxes without outline), sub-steps of the methodology within the main steps (light grey boxes with black dashed outline), data collected (white boxes with black dashed outline), methodology results (white boxes with red solid outline).

Running the InVEST Model
The data required to run the InVEST SWYM (Figure 1 and Appendix A, Table A1) were collected for a 15-year time period, spanning 1 January 2000 to 31 December 2014.Data were gathered to allow for the analysis of watersheds ranging in size from 174 km 2 to 10,814 km 2 with a geographic distribution covering most of North America between 26.9 • and 60.9 • N and 53.6 • and 137.1 • W. We used the daily Climate Prediction Center (CPC) Global Unified Precipitation and Temperature data (https://www.esrl.noaa.gov/psd/)and MOD16 monthly reference evapotranspiration rasters from the Numerical Terradynamic Simulation Group (NTSD; http://www.ntsg.umt.edu) for our analyses.The climate data were the only data that varied with each run of the model.All other data required by the InVEST SWYM stayed constant for the duration of the runs.This allowed us to use single global rasters for soils [22], the digital elevation model (DEM), and land use/land cover [23].The soil data were available at seven different depths ranging from 0 cm to 200 cm.For our analyses, we selected the 30 cm depth as we felt it was a reasonable transition point between runoff versus groundwater recharge.A quick comparison between the surface layers (0 cm and 30 cm) found that, in general, they had similar soil classifications.For the DEM data, we used the HydroSHEDS [24] hydrologically conditioned elevation data.The individual 5 × 5 tiles were stitched together using GDAL Warp [25] to produce a single North American DEM.See Appendix A, Table A1 for a complete list of data sources and attributions.
We selected watersheds for inclusion in our analysis based on a list of 1545 gauging stations with validated flow data (available from The Global Runoff Data Centre, 56,068 Koblenz, Germany (GRDC)).Gauging station data were limited to stations in North America with daily flow data from 1 January 2000 to 31 December 2014.The GRDC also provides watershed boundaries for each of the gauging stations.That is, each gauging station is the pour point for the associated watershed shapefile [26].We used each gauging station specific watershed shapefile to clip the global raster files (e.g., DEM, precipitation, temperature, soils, evapotranspiration and land use/land cover) during each run of the InVEST model to ensure that the results were specific to the watershed being modeled.The resulting 749 modeled watersheds represent a broad range of environmental conditions (e.g., annual precipitation, the annual number of days with below freezing temperatures, thermal amplitude, and altitudinal difference), which can be observed in Appendix A, Table A1.
The InVEST SWYM model requires monthly total precipitation and not daily totals; however, it also requires the number of precipitation events per month for the watershed that is being modeled as data input.To calculate the number of precipitation events per month, we averaged the total daily precipitation across all pixels contained within the boundaries of each GRDC watershed shape file.Days with an average precipitation per pixel of greater than 0.1 mm were counted as precipitation events and monthly totals were tabulated for inclusion in the SWYM run.Due to the requirement for monthly precipitation, we converted daily total precipitation to monthly total precipitation by summing the daily precipitation values at each pixel of the raster to create monthly precipitation rasters.
Running the InVEST SWYM model for a single watershed or set of parameters is easily accomplished through the graphical user interface of the application; however, as we wanted to run the model for 749 watersheds over a 15-year period we used the InVEST Python application programming interface (API) (Figure 1).Details of the API are available from http://invest.readthedocs.io/en/latest/.Using the InVEST API allowed us to easily clip the rasters for each watershed to the boundaries and to re-project it to the relevant UTM based coordinate reference system (CRS) of each watershed prior to running the model.It also allowed us to programmatically update the parameters so that we could step through all combinations of watershed × year.As part of the re-projection step, we also resampled the DEM rasters to a pixel size of 90 × 90 m because the InVEST model expects the pixels of the DEM raster to have sides of equal length.The InVEST SWYM has several parameters that can be optimized to improve model performance (α, β, γ and flow accumulation threshold).As we ran the model for 749 watersheds, it was not feasible to optimize these parameters for each watershed.As such we used the default settings from the SWYM for α, β and γ (α = 1/12, β = 1.0, γ = 1.0).We did set the threshold value to 125, which coincides with a 1 km 2 contribution area.Explanations of each of the parameters can be found on the Natural Capital Projects website (http://data.naturalcapitalproject.org/nightlybuild/invest-users-guide/html/seasonal_water_yield.html).

InVEST Model Output
Among several other model outputs, the InVEST SWYM produces 12-monthly raster files (per year run) of predicted total monthly Quick Flow (QF or runoff) per pixel (Figure 1).The total predicted monthly runoff (QF) for the watershed is then calculated as the sum of all pixels within the monthly QF raster.We extracted these monthly predicted QF values for every combination of watershed × year × month.This provided us with a set of monthly predicted flow for each watershed over the 15-year period.These data were then compared to GRDC gauging station data of observed monthly flow for each watershed (Figure 1).

Comparative Analysis of Predicted versus Observed Flow
To compare the predicted flow from the InVEST SWYM model to the observed flow of the gauging station data, we squared the Pearson's product-moment correlation (r) for each watershed (Figure 1) to create a simple standardized measure of model performance.A squared Pearson's product-moment correlation allowed us to focus on the strength and not the direction (positive or negative) of the predicted flow patterns.As we were not examining the significance of r, normality was not required in the underlying data.This straightforward approach allowed us to determine whether the predicted flow from the SWYM co-varied with the actually observed flow for each month of the study period for each watershed and to test for regional patterns in SWYM performance with respect to broad-scale climate variables.

Defining Regions and Characterization of Climate Variables for Each Region
We used a Principal Component Analysis (PCA) and a Hierarchical Cluster Analysis (HCA) to examine how the model results were influenced by climate (minimum and maximum temperature, precipitation, vapor pressure, solar radiation, and days below zero) and geographical variables (mean, maximum, and minimum elevation) for each watershed (Figure 1).This combined analytical approach has been shown to be efficient in uncovering patterns in climatic and geographic variables (e.g., References [27,28]).Before running the PCA, we standardized all environmental variables of the watersheds (Appendix B, Table A2) using the "scale" function in base R. We did this to ensure that all the terms were non-dimensional and the significant differences in the magnitude between variables were minimized.We then used the PCA analysis from the "vegan" package [29] to summarize the significant environmental variability.We used significant PCA axes selected by the Broken Stick method [30], using the "bstick" and "screeplot" functions of "vegan" [29].Subsequently, we retained only these significant axes of the PCA site scores to perform the HCA based on Euclidean distances.For this, we employed the "agnes" function of the "cluster" package [31].We visually assessed the dendrogram to determine which tree produced distinct groups and minimized noise when plotting.We used the "cutree" function of the "cluster" package [31] to produce the groups.All analyses were performed using R version 3.4.0[32].We created boxplots of the environmental variables for each bioclimatic region, which allowed us to explore the variation between regions and to visualize five summary statistics (the median, two hinges (first and third quartiles) and two whiskers).We also created boxplots for the results for each region of the InVEST SWYM.
Finally, we performed an AIC-based multiple regression to analyze which environmental variables (independent variables) affected the results of the hydrological model in each bioclimatic region (measured above as r 2 between observed and predicted flow; dependent variable) (Figure 1).The AIC analysis allowed us to determine the most parsimonious group of variables that best predicted the fit of the model for each environmental region.Once defined, we conducted a targeted multiple regression analysis to identify which of the AIC-derived variables explained most of the variance in the model fit (Figure 1).Each independent variable in the multiple regressions had an associated beta coefficient (β) indicating how much the standard deviation of the dependent variable increased when the independent variable was increased by one standard deviation (SD), assuming other variables in the model were unchanged.The net values of β are a measure of the total effect of the independent variables, so the independent variable with the highest net value of β is the one with the greatest total effect over the dependent variable.Thus, in this study, the analysis of net values of β showed the relative importance of the environmental variables for the variability of r 2 between observed and model-predicted flow.We also used partial correlations to verify our results.This analysis was carried out using the software Statistica 7.

Improving Model Performance
After analyzing the InVEST SWYM results and determining the key environmental factors associated with each of the bioclimatic regions, we then carried out a case study (Figure 1) to re-evaluate model performance when additional parameters were included.To do this, we selected the bioclimatic region containing the most watersheds in this study.This region, the Low elevation Boreal Sub-humid (LeHSh) one, was also the region in which the InVEST SWYM had the lowest mean r 2 values (Figure 1).The environmental variables associated with this region suggested cryospheric conditions might be an important factor.
The InVEST SWYM assumes that precipitation falls as rain and is immediately available for partitioning into the different hydrological pathways [15].However, precipitation that falls as snow or ice can accumulate and act as a reservoir in watersheds that experience periods of below freezing temperatures [33].The accumulated ice and snow can be released either gradually, as a single pulse, or be locked away for prolonged periods of time (e.g., glaciers).Consequently, we coupled a snow accumulation and ablation model to the InVEST Seasonal Water Yield Model (SWYM + SNOW17) to account for the effect of cryospheric conditions when using the water yield model [34].This resulted in a second set of predicted flows for each watershed in the LeHSh region that had been cryospherically conditioned by the SNOW-17 model.The predicted flow patterns from the SWYM + SNOW-17 model were compared to the observed flow patterns of the gauging station data using the squared Pearson's product-moment correlations for each watershed (Figure 1).
The SNOW-17 snow accumulation and ablation model [34] uses temperature and precipitation data to model how much precipitation accumulates as snow or ice versus falling as rain.It also models the conversion of accumulated snow or ice to either water or water vapor.For the purposes of our study, the SNOW-17 model was used to determine how much precipitation fell as rain, which was available to the InVEST SWYM immediately; how much accumulated as snow and was unavailable to the InVEST SWYM; the conversion of accumulated snow to water, which was available to the InVEST SWYM for flow partitioning versus being lost to evaporative processes.The additional data requirements of the SNOW-17 model were easily satisfied by using the daily climate data already collected for the InVEST model as discussed in the analysis of the fit of the hydrological model section above.

Analysis of the Fit of the InVEST SWYM
When analyzing 749 watersheds from across North America, we found the correlations between modeled flows and actual flows to be highly dispersed and relatively poor.Ninety-two percent of the r 2 values were less than 0.5 and 42% were less than 0.1 (Figure 2).

Analysis of the Fit of InVEST SWYM by Region
The first two components of the PCA analysis explained 81% of the total variance within the watersheds (Figure 3 and Table 1).The first PCA axis explains 53% of the variance (Figure 3C and Table 1) and is associated with climatic variability; positive scores occur for vapor pressure, minimum, and maximum temperature, whereas negative scores were associated with the number of days with below zero temperatures (Figure 3C and Table 2).The PCA axis 2 explains 28% of the variance and is mainly associated with elevation (Figure 3B and Table 1); positive loadings occur for mean, minimum, and maximum elevation (Figure 3C and Table 2).Thus, we interpret axis 1 of the PCA as representative of climatic variability and axis 2 as representative of a topographic gradient.
The cluster analysis based on the two first components of the PCA allowed us to identify five distinct environmental regions across our study area (Figure 3A,B).Their geographic distribution can be observed in Figure 4 and their environmental characteristic in Figure 5. Clusters 1 and 2 comprised all watersheds under Boreal Sub-humid conditions.Cluster 1 grouped 108 watersheds on the North Western mountain areas of the continent (Mountain Boreal Sub-humid, Figures 3A and 4     The InVEST SWYM model produced better results in some regions (LeTH and MeTHh) than in others (MBSh, LeBSh, and MTA) (Figure 6).The r 2 value between actual and modeled flow in LeTH was 0.3 ± 0.01 (mean ± standard error, respectively), while in the MeTHh the r 2 was 0.35 ± 0.05 (Figure 6).These watersheds are under temperate, humid or hyper-humid climate conditions and mid to low elevation lands (Figure 5).The r 2 values between actual and modeled flow were low in LeBSh (0.16 ± 0.01), MBSh (0.17 ± 0.2), and MTA (0.16 ± 0.2) (Figure 6).In contrast to the former two environmental regions, where the hydrological model produced better results, the LeBSh and MBSh were under cold climate conditions, while the LeBSh and MTA were mountainous regions (Figure 5).In addition, the r 2 values between observed and modeled flow were significantly correlated with different environmental variables in each of the five climatic regions (Table 3).

Analyzing Environmental Features to Improve InVEST SWYM Results
The watersheds in the LeBSh environmental region had the lowest mean r 2 values of the correlation between actual and modeled flows (Figure 6).This region contained climatic features related to cryospheric conditions (Figure 5).The watersheds in this region had more than 120 days with below zero temperatures and the second highest number of days with below zero temperatures between the five environmental areas defined in this study, with temperatures reaching −20 • C (the lowest compared with the other four environmental groups) (Figure 5).
For LeBSh, the AIC-based multiple regression indicated that ∼12% (adjusted r 2 = 0.12; p < 0.001) of the variance of the model fit (measured as r 2 between observed and model-predicted flow) could be accounted for by a linear combination of the independent variables defined by the AIC analysis (minimum and maximum temperature, precipitation, days below zero, maximum and minimum elevation, and vapor pressure) (Table 3).Net values of the beta coefficient (β) show the importance of the variables in a linear multiple regression.The variables that best explained the model fit for LeBSh were the number of days with below zero temperatures (β = −0.18)and the minimum temperature (β = −0.57);both related by cryospheric conditions (Table 3).Partial correlation shows the same results as the net values of the β (Table 3).Variables related to topographic features were also important in explaining the model fit in this region (minimum elevation, β = 0.53; maximum elevation, β = 0.46).The number of days with below zero temperatures was negatively correlated with the model results.
This indicates the model performance degrades as the number of days of below freezing temperatures increases (Table 3).In addition, the minimum temperature was negatively correlated with the model results, indicating that when minimum temperatures increase, the model results were poorer.
When comparing the LeBSh with the MeTHh group (which is the environmental region where the model produced the best results), we found that the LeBSh region had poorer model predictions (r 2 = 0.16) than the MeTHh group (r 2 = 0.35) (Figure 6).For MeTHh, the AIC-based multiple regression indicated that ∼80% (adjusted R 2 = 0.8; p < 0.001) of the variance of the model fit could be accounted for by a linear combination of the independent variables defined by the AIC analysis (precipitation, maximum temperature, solar radiation, and vapor pressure) (Table 3).The variable that best described the model fit for MeTHh was the maximum temperature (β = 0.67) (Table 3).The maximum temperature was positively correlated with the model results, indicating that when maximum temperatures increased, the model results improved.
In addition, in the MeTHh region, the mean number of days with below zero temperatures was 53 (67 days less than in LeBSh) and the mean minimum temperature was −5.5 • C (14.5 • C higher than in LeBSh) (Figure 5).The InVEST SWYM, therefore, performed better in regions with a fewer number of days of below zero temperature and a higher minimum temperature.These results suggest that cryospheric conditions might be influencing the results of the hydrological model.
After accounting for cryospheric conditions and running the SWYM (SWYM + SNOW-17) for the watersheds located inside the LeBSh region, correlation values between the actual and modeled flows improved in 30% (114 watersheds) of the watersheds (Table 4).These improvements can also be observed in a shift of the histogram distribution, with considerably fewer watersheds in the first bin (0 to 0.1) and increased in all other bins (Figure 7).The mean r 2 value between real and modeled flow increased by 10%.

Discussion
Our assessment of the performance of the InVEST SWYM across the differing environmental regions demonstrated that the correlation between the modeled and actual flows was highly variable and relatively poor.To investigate the regional differences further, we grouped the watersheds into five different environmental regions within North America, similar to those proposed for North America by Rivas-Martinez [35] and evaluated the differences in model performance to identify key variables.The results of this analysis indicated that identifying features such as bioclimatic regions and accounting for these differences can substantially improve the model output.We readily demonstrated this by examining performance gains due to the incorporation of cryospheric conditions associated with the region with the poorest model performance (the Low elevation Boreal Sub-humid region-LeBSh).In this region, the SWYM results were significantly improved by the addition of cryospheric variables.
After we applied a snow accumulation and ablation component (SNOW-17) to the InVEST SWYM, runoff predictions significantly improved in 30% of the watersheds in the LeBSh region.Our results, thus, highlight that the benefits of incorporating additional processes, such as the SNOW-17 accumulation and ablation model, into the simple hydrological model can be large relative to the cost of implementing them.For example, the SNOW-17 model required very little additional data in the form of daily versus monthly temperature and precipitation data to run.These data are available for most regions of the world from publicly available sources (see Appendix A, Table A1).Although the SNOW-17 model did not require a great deal of additional data, it does take into account most of the physical processes that occur within snow cover [34] and implementing the model was straightforward with available existing code sources available (e.g., https://github.com/UW-Hydro/tonic/blob/master/tonic/models/snow17/snow17.py).
It is worth noting that the correlation between the predicted flows using the InVEST SWYM and the actual flows was overall relatively poor, with 90% of the r 2 values less than 0.5.In regions under temperate climate conditions, the correlation between the predicted and actual flows was less than 0.5 in 77% of the watersheds.Clearly, variables other than those related to cryospheric conditions were affecting hydrological model results.As an example, in our study, vapor pressure and variables related with terrain slope (e.g., mean, min and max elevation) were highly correlated with the hydrological results in some specific environmental regions (MBSh, and LeBSh).In other studies, landscape characteristics related to terrain topography (e.g., slope) have been shown to produce major changes in the hydrological model results [36].Furthermore, for three environmental regions (LeBSh, LeTh, and MeTHh), the environmental variables used in this study explain a low percentage of the variance in model fit suggesting that other variables (e.g., land cover, dams, soil characteristics, etc.) might be affecting hydrological model results.Under such conditions, the model may need to be expanded or to consider specific parameters that were not originally incorporated here.
Our results are similar to other studies that highlight the importance of climatic inputs for hydrological models [37,38].In particular, the hydrological model we used (InVEST SWYM) is most sensitive to uncertainty in climate forcing such as precipitation [39][40][41].As described by Hamel and Guswa [41], errors in climate data are spatially heterogeneous and can lead to differences in the model results in diverse environmental regions, as was the case here.In that regard, there are other global products that provide monthly climate information in raster format (e.g., CHELSA; http://chelsa-climate.org/).It would be worth comparing the results from the InVEST SWYM using another input dataset to investigate if poor model results are related to data source issues.It is clear from our work, however, that accounting for unaddressed hydrological processes, or by improving existing processes in the InVEST SWYM, significant gains in model performance can be made.
Optimizing available model parameters for each watershed may help to account for the variation within a watershed and improve model performance; however, as mentioned previously, this can be a time consuming and difficult process.In situations where reliable data do not exist, model calibration may not be possible.As such, improving the base model by identifying critical limitations through this type of analysis (versus optimizing for individual watersheds) may prove to be more broadly beneficial and applicable.We have assumed that the global datasets that we have used (precipitation, soils, and flow data) are all reasonably accurate for each of the watersheds included in the study and acknowledge that they may not be the best data to use; however, as mentioned above, some regions are data poor and these are likely to be the data sets that they turn to for such an analysis.Furthermore, improving the models based on broadly available data is likely to benefit more end users than to focus on improvements at the level of the individual watershed.
While our results seem to indicate that the InVEST SWYM performs poorly, it should be made clear that our analyses were not designed to address this question.We did not explicitly test the performance of the SWYM.That is, we did not optimize the run parameters (α, β, γ and flow accumulation threshold) for any of the watersheds and used the model with all of its default settings, as optimizing for each watershed was not feasible.More importantly, we were looking for indicators of broad-scale environmental factors that influence hydrologic processes that the model may not incorporate.Our multiple regression results (Table 3) give some indication of what these variables might be.Days below zero clearly indicated the importance of cryospheric conditions but the consistent significant values of variables related to evaporation processes (such as vapor pressure and temperature) in the regression results indicated that evapotranspiration in certain environmental regions might be inadequately addressed by the InVEST SWYM (Table 3).To make a robust assessment of the InVEST SWYMs performance, a careful benchmarking study that involves watershed-specific variable optimization is required.
Finally, it should be acknowledged that we made several efforts to include Central and South America's watersheds in this study with little success.We limited our data acquisition to watersheds with both flow and climate data covering the 15-year period between 1 January 2000 and 31 December 2014.This excluded all but two watersheds for South America (SA) as most data for this region did not comply with our criteria.When analyzing the two South American watersheds, we found they had environmental conditions similar to those of the MTA region (northern watersheds in Figure 8) and the LeTH region (southern watersheds in Figure 8).The model fits (r 2 ) in these watersheds were 0.33 and 0.24 respectively.However, due to the lack of SA sites, it was not possible to specifically discuss how the InVEST model performed in this part of the continent but many of the same hydrological and climate variables may be important in this region where conditions are similar.However, as mentioned previously, hydrological models have uncertainty related to model structure and their performance might not be appropriate if the model structure is based on incorrect dominant processes.Thus, it is necessary to test the transferability of the model structures under the specified environmental conditions in a particular region.When comparing the Northern with the Southern Hemisphere, the former is water-rich and contains the majority of the world's fresh water with 60% of watersheds experiencing cryospheric conditions versus the Southern Hemisphere, which is water poor, tends to be more arid and has significantly fewer lakes and rivers [42].Including South America's watersheds could provide greater insights into important hydrological processes by providing additional data characterizing conditions that are not overwhelmingly moist and cool.However, as demonstrated by our results, running this model with global data is not always possible in SA, further highlighting the problem of obtaining quality data for these types of analyses and the need for improved data availability in data-poor regions of the world.

Conclusions
Our analyses demonstrated that the correlation between actual flow patterns and predicted flow patterns from the Natural Capital Project's InVEST SWYM varied depending on the environmental conditions.However, by isolating the environmental factors that caused the model to generate poorer results, we were able to identify important environmental variables and improve the performance of the model.For example, in the Low elevation Boreal Sub-humid region (LeBSh) of North America, the model produced the poorest results and these were correlated with the cryospheric conditions in the region (days below freezing temperature and minimum temperature).After we applied a specific cryospheric model (SNOW-17) to the SWYM, runoff predictions significantly improved in the 30% of the watersheds included in the aforementioned region, underlining the need to incorporate this and other standard models of watershed processes into the InVEST SWYM framework whenever possible.
We also highlight the importance of understanding watershed-specific climatic conditions before implementing a hydrological model.Optimization and calibration should be undertaken whenever possible, even if data is only available for another watershed in the same region, to address variables that may be influential for those particular conditions.We have shown that accounting for cryospheric conditions improved hydrological model results in an environmental region under low elevation, boreal sub-humid climate.However, other important climatic conditions, such as evapotranspiration, solar radiation, and slope are usually oversimplified in some hydrological models and would need specific attention.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 2 .
Figure 2. Histogram of the watersheds grouped by the r 2 values of the correlation between actual and modeled flow.The interval size is 10%.
) and cluster 2 grouped 383 basins on low lands of the North Central and North Eastern regions of North America (Low elevation Boreal Sub-humid, Figures 3A, 4and 5).Clusters 3, 4, and 5 comprised all watersheds under temperate conditions.Cluster 3 is composed of 165 basins located in low elevation humid areas in Southeastern North America (Low elevation Temperate Humid, Figures 3A, 4 and 5).Cluster 4 grouped 73 basins on mountain arid areas on the Center and Southwest of the continent (Mountain Temperate Arid, Figures 3A, 4 and 5).Finally, cluster 5 included 22 watersheds under Oceanic Hyper-humid conditions on mid-elevation areas in the Northwest of the continent (Mid-elevation Temperate Hyper-humid, Figures 3A, 4 and 5).

Figure 4 .
Figure 4. Geographic distribution of watersheds included in the study.The color of the dot indicates the environmental region each watershed was grouped into by the cluster analysis based on environmental conditions.

Figure 6 .
Figure 6.Mean (Inner Square), standard error (bar), 95% confident interval (whisker) of the r 2 values of the correlation between actual and modeled flow in each of the environmental regions.

Figure 7 .
Figure 7. Histograms of the correlation values (r 2 ) between actual and modeled flow with the SWYM (A) and with the SWYM+SNOW17 (B).(C) Mean (Inner Square), standard error (bar), 95% confident interval (whisker) of the correlation between actual and modeled flow in the LeBSh environmental regions with SWYM and SWYM + SNOW17.

Figure 8 .
Figure 8. Geographic localization of the two South American watersheds included in the study.

Table 1 .
Importance of components.

Table 2 .
Factor loadings of the environmental variables on Principal Component (PC) axes 1 (PC1) and 2 (PC2), expressed as correlations between a variable and a particular axis.

Table 3 .
AIC-based multiple regression analysis of the r 2 between observed and model-predicted flow in the environmental regions.Adjusted r 2 of the multiple regression, β coefficient, and partial correlations (in brackets) for the variables used in the multiple regression.β coefficients and partial correlations values marked in red are significant at p < 0.05.

Table 4 .
Number and percentages of watersheds where the correlation between the actual and modeled flows improved (Fisher r to Z transform, α = 0.05) after running the SWYM accounting for cryospheric conditions (SWYM + SNOW-17).

Table A1 .
Data types and sources required to run the InVEST model.

Table A2 .
Raw data used for PCA analysis.