Evaluation of InVEST’s Water Ecosystem Service Models in a Brazilian Subtropical Basin

Abstract

The biophysical modeling of water ecosystem services is crucial to understanding their availability, vulnerabilities, and fluxes. Among the most popular models, the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) models stand out. While many studies have used them, few have assessed their performance. This study evaluates the performance of InVEST’s Seasonal Water Yield, Nutrient Delivery Ratio, and Sediment Delivery Ratio models in a subtropical basin in southeastern Brazil on temporal and spatial scales, using 39 years of streamflow data, 29 for total phosphorus and total nitrogen, and 19 for total suspended solids. Statistical indicators R², PBIAS, and NSE, were also calculated. The performance of the models varied according to the type of simulated WES and analysis scales used, with the Seasonal Water Yield model demonstrating the best performance and effectively representing the spatial and temporal variability of the average annual streamflow. All models performed well in simulating long-term mean values when compared to observed data. While one should bear in mind the study’s limitations, the results indicate that the models perform well in terms of relative magnitude, although their application in studies involving water-resource management and decision making is limited.

1. Introduction

A water ecosystem services (WES)-based approach, incorporating an integrated analysis among multiple pressures, ecological statuses, and ecosystem services has been identified as an effective tool to plan and manage water resources, as it links environmental conservation with socioeconomic development [1]. Such an approach enables an assessment of how anthropic activities affect ecosystem composition and functioning, impacting WES, including water supply services for human use, irrigation, and energy generation, as well as regulation services associated with flood mitigation, erosion control, and water purification [2].

A pillar of this approach is biophysical modeling, which quantifies WES from mathematical equations in computational environments [3]. Biophysical models facilitate the understanding of ecohydrological processes in hydrographic basins and the assessment and forecasting of multiple pressures [1]. Well-calibrated, validated models have numerous applications in water management, such as quantifying the effects of land use and climate change [4,5,6], assessing flood risk [7], and designing and evaluating forest restoration programs [8,9,10,11], in addition to generating and transferring historical data for unmonitored basins, which is critical for hydrological studies [12]. Such applications can help formulate and evaluate environmental conservation policies that support sound decision making in water resource management and land-use planning [13].

Over the past few decades, free models to map ecosystem services have emerged [14,15,16], ranging from simple approaches based on land-use data or habitat-based proxies to complex models that quantify physical processes in ecosystems [3]. Traditional hydrological models such as the Soil and Water Assessment Tool (SWAT) [17], Variable Infiltration Capacity model [18], and the Regional Hydro-Ecological Simulation System [19] can be used to map and quantify WES. While they can accurately represent ecohydrological processes critical to a broad spectrum of ecosystem services, such as drinking water and recreation [1,20], their application is limited in regions with a dearth of data since their extensive parameters require detailed data about the study area as well as users with expertise in hydrology and modeling for calibration and validation [12,21].

Among simpler alternative models such as Artificial Intelligence for Ecosystem Services [14], Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) [16], and the Multiscale Integrated Model of Ecosystem Services [15], the InVEST’s models stand out and are frequently used to quantify and map WES. InVEST works with hydrological simplifications at the watershed scale that enable the simulation of such services as annual and seasonal water yield, water purification, and erosion control [16]. InVEST’s models have been applied in diverse WES contexts and locations, including assessing the effects of mitigation policies in Europe [22], climate change and land use in WES in the United States [5], and ecological restoration in China [23]. Although these models are applied globally, there are few studies that investigate their performance in comparison with empirical values [13,24,25,26,27,28,29]. Of these, only Lu et al. [29] analyzed the temporal performance of the Seasonal Water Yield (SWY) model, while the vast majority evaluated the spatial performance of InVEST models using long-term average values.

Assessing the performance of ecohydrological models is crucial to understanding the significance of parameters and input data to the reliability of their results, which inform decision-making and policy formulation [30]. In addition, such assessments can aid developers seeking to enhance the effectiveness and efficiency of such tools [26] and users looking to understand their potential and limitations in watersheds [31], since their performance may vary according to regional climatic and hydrogeological characteristics.

A limited number of studies address the use of InVEST’s WES models in Brazil [32,33,34,35,36,37,38], but none analyze their performances, comparing them with observed data. This study sought to calibrate and validate InVEST’s SWY, Nutrient Delivery Ratio (NDR), and Sediment Delivery Ratio (SDR) models, considering their ability to simulate WES spatial and temporal variability. The study investigates the hypothesis that InVEST’s WES models can effectively represent the observed annual values of streamflow, exported sediment, total nitrogen (TN), and exported total phosphorus (TP) in a subtropical watershed located in southeastern Brazil. The specific objectives were to assess the models’ sensitivity to variation in calibration parameters, calibrate and validate the models, and evaluate the performance of the models. It is hoped that the results of this study can clarify the strengths and limitations of InVEST WES models to assist those using these tools to formulate policies for the conservation of water resources.

2. Materials and Methods

2.1. Study Area

The study was carried out in the Jacaré-Guaçu River Basin, in the central-eastern portion of the state of São Paulo, between the geographical coordinates 21°37′00″ S and 22°22′00″ S and 47°43′00″ W and 48°56′00″ W (Figure 1). The drainage area is nearly 4172.12 km², covering the municipalities of Araraquara, Boa Esperança do Sul, Gavião Peixoto, Ibaté, Ibitinga, Itirapina, Nova Europa, Ribeirão Bonito, São Carlos, and Tabatinga, and is a sub-basin of the water resources management unit of the State of São Paulo (UGRHI 13), viz., the Tietê-Jacaré Hydrographic Basin. The Jacaré-Guaçu River is a direct tributary of the Tietê River, the largest watercourse in the state of São Paulo, and a source of water for UGRHI 13. It is formed by the junction of the Lobo and Feijão streams, which originate in the Itaqueri mountain range and Cuscuzeiro, covering 155 km until reaching the Ibitinga reservoir in the mid-course of the Tietê River [39]. The basin includes small reservoirs that produce energy, such as the Lobo and Santana plants, and others that supply small rural properties.

The Jacaré-Guaçu River Basin lies between the CWA and CWB climatic zones, with dry winters and humid summers [40]. Average annual precipitation is approximately 1400 mm, with a rainy season between October and March (Figure 2), which accounts for an average 80% of the total annual precipitation. Average monthly temperatures in the region range from 17.4 °C in June to 23.7 °C in February.

The topographic elevation varies from 372.0 m to 1024.0 m above sea level. The basin has a predominantly flat relief, composed of smoothly undulating hills (around 80%), arising from weathering of the Botucatu and Pirambóia formations, with low drainage density. The high slopes are mainly concentrated in areas of undulating and steep relief, near the sources of the Jacaré-Guaçu River and the drainages in the central portion of the basin.

The basin is on the Guarani Aquifer System (GAS), formed by Aeolian sandstones from the Jurassic (Botucatu Formation) and fluvio-eolic Triassic (Pirambóia Formation) periods, which form underground water reservoirs. GAS is one of the world’s largest underground freshwater reservoirs, with 1,200,000 km² covering Brazil, Argentina, Paraguay, and Uruguay [41]. Aquifer recharge occurs in approximately 10% of the upwelling areas of the Botucatu and Pirambóia formations [42]. The basin is an important GAS recharge area, where the aquifer is unconfined. At the basin, the aquifer covers an extensive area (~1640 km²) ranging from the region of Itirapina, at the head of the basin, to Ibitinga, close to its outlet. In addition to GAS, the Bauru and Serra Geral aquifers are found in the basin [39].

The basin features 10 types of soils: Red-Yellow Argisol, Haplic Gleysol, Red-Yellow Latosol, Red Latosol, Haplic Luvisol, Litholic Neosol, Quartzarenic Neosol, Red Nitosol, Haplic Organosol, and Haplic Planosol. Latosols constitute the predominant class and cover 60% of the basin, followed by Red-Yellow Argisols (14.77%) and Quartzarenic Neosols (12.77%) [43].

The basin’s native vegetation is composed of grassland and forest, represented by remnants of Brazilian Cerrado (Savanna), riparian, and semi-deciduous forests [39]. Intense changes have been observed in the landscape as a result of expanded agricultural activities, in particular, sugarcane cultivation, that have significantly reduced natural vegetation [44]. At present, natural areas occupy only 17% of the basin [45]. The reduction of the basin’s natural areas, such as wetlands and riparian forests adversely impacts the region’s biodiversity and ecosystem services.

2.2. Databases

The SWY, NDR, and SDR models were used with data from monitoring stations and high spatial resolution cartographic data. ALOS-PALSAR images with a spatial resolution of 12.5 m, freely available at the Alaska Satellite Facility [46], were used to determine DEM. Land-use data from 1985 to 2019 with a spatial resolution of 30 m were obtained from the Brazilian Annual Land Use and Land Cover Mapping Project (5th collection) [45]. Meteorological data were attained directly from the climatological station of the Center for Water Resources and Environmental Studies of the School of Engineering of São Carlos. Precipitation data from the monitoring stations of the National Water and Sanitation Agency (ANA) [47] and data estimated by the Climate Hazards Group Infrared Precipitation with Stations of high resolution (0.05°) obtained via Google Earth Engine were also used. Daily and monthly streamflow data were obtained from ANA and the Department of Water and Energy of the State of São Paulo (DAEE) [47,48]. Data on total solids (TS), dissolved solids (DS), and TN and TP concentrations were acquired from the Environmental Company of the State of São Paulo [49]. Data characterizing the basin’s soils were extracted from the pedological map of the state of São Paulo, at a scale of 1:100,000 [43].

2.3. Water Ecosystem Service Models

InVEST’s ecosystem services models are free, open-source software used to evaluate a broad range of ecosystem services [16]. They are spatially explicit, that is, they use maps as input data and produce them as output data. The operating structures of each model used in this study are detailed below. Further information can be found in Sharp et al. [16]. The biophysical tables and all parameters required to run InVEST’s WES models are available in the supplementary biophysical tables (Tables S1 and S2).

2.3.1. Seasonal Water Yield Model

The SWY model was used to quantify annual water production in the basin. The algorithm enables quantifying annual base flow (BF) and annual and monthly surface runoff (QF) for each pixel. The sum of BF and QF characterizes the basin’s annual streamflow. Monthly QF is calculated based on a modified approach of the Natural Resource Conservation Service (NRCS) curve number method [50], which determines it from monthly rainfall and the number of rainfall events (Equation (1)).

$${\mathrm{QF}}_{\mathrm{i},\mathrm{m}}={\mathrm{n}}_{\mathrm{m} }\left(\left({\mathrm{a}}_{\mathrm{i},\mathrm{m}}-{ \mathrm{S}}_{\mathrm{i}}\right)\exp \left(-\frac{0.2{\mathrm{S}}_{\mathrm{i}}}{{\mathrm{a}}_{\mathrm{i},\mathrm{m}}}\right)+\frac{{\mathrm{S}}_{\mathrm{i}}^2}{{\mathrm{a}}_{\mathrm{i},\mathrm{m}}}\exp \left(\frac{0.8{\mathrm{S}}_{\mathrm{i}}}{{\mathrm{a}}_{\mathrm{i},\mathrm{m}}}\right){\mathrm{E}}_{\mathrm{i}}\left(\frac{{\mathrm{S}}_{\mathrm{i}}}{{\mathrm{a}}_{\mathrm{i},\mathrm{m}}}\right)\right)\times \left(25.4 \left[\frac{\mathrm{mm}}{\mathrm{in}}\right]\right)$$

(1)

$${\mathrm{a}}_{\mathrm{i},\mathrm{m}}=\frac{{\mathrm{P}}_{\mathrm{i},\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}}/25.4$$

(2)

where P_i,m is the monthly precipitation in pixel i for month m (mm), n_m is the number of rain events in pixel i in month m, a_i,m is the mean rain depth on a rainy day at pixel i in month m (Equation (2)), S_i = 1000/CN_i − 10; CN_i is the curve number (Table S1) for pixel I, and E_i is the exponential integral function. The value 25.4 is an inch-to-millimeter conversion parameter. The annual QF_i is the sum of the monthly values QF_i,m (mm).

The local recharge for each pixel is calculated from the local water balance (Equation (3)). The local recharge value is determined on an annual scale, using monthly values. The portion of available monthly rainfall that is not runoff over the surface is divided between local recharge and evapotranspiration. For a given pixel, the partitioning is affected by upgradient recharge, and parameters α and β influence the availability of groundwater for evapotranspiration (Equation (4)).

$${\mathrm{L}}_{\mathrm{i}}={\mathrm{P}}_{\mathrm{i}}-{\mathrm{QF}}_{\mathrm{i}}-{\mathrm{AET}}_{\mathrm{i}}$$

(3)

$${\mathrm{AET}}_{\mathrm{i},\mathrm{m}}=\min ({\mathrm{PET}}_{\mathrm{i},\mathrm{m}}; {\mathrm{P}}_{\mathrm{i},\mathrm{m}}-{\mathrm{QF}}_{\mathrm{i},\mathrm{m}}+{ \mathsf{\alpha }}_{\mathrm{m}}{\mathsf{\beta }\mathrm{L}}_{\mathrm{Sum}.\mathrm{avail};\mathrm{i}})$$

(4)

$${\mathrm{PET}}_{\mathrm{i},\mathrm{m}}={ \mathrm{K}}_{\mathrm{c},\mathrm{i},\mathrm{m}} .{ \mathrm{ET}}_{0,\mathrm{i},\mathrm{m}}$$

(5)

where L_i is the local recharge for pixel i (mm), P_i is the annual precipitation for pixel i (mm), AET_i is the real annual evapotranspiration in pixel i (mm), AET_i,m is the real monthly evapotranspiration in pixel i (mm), PET_i,m is the monthly potential evapotranspiration at pixel i (mm), αm is the fraction of upslope annual recharge available in month m (default: 1/12), β_i is the fraction of the upslope subsidy available for downslope evapotranspiration (default: 1), L_Sum.avail,i is the sum of upslope subsurface water potentially available at pixel i, K_c,i,m is the monthly crop coefficient for ground cover type (Table S1) in pixel i, and ET_0,i,m is the reference evapotranspiration (Equation (S1)) at pixel i for month m (mm).

The annual base flow index (B), which represents the portion of water that reaches the watercourse via groundwater, is calculated based on local recharge values. Negative values indicate that the pixel does not contribute to BF, and B is assigned a value of zero. When the pixel contributes to BF, the value of B is calculated as a function of the amount of flux leaving the pixel and its relative contribution to BF reloading this pixel. For a pixel that is not adjacent to the stream channel, the cumulative base flow (B_sum,i) is proportional to the cumulative base flow that is leaving the adjacent downslope pixels minus the cumulative base flow that was generated on that same downslope pixel.

$${\mathrm{B}}_{\mathrm{sum},\mathrm{i}}={ \mathrm{L}}_{\mathrm{sum},\mathrm{i}}, \mathrm{if} \mathrm{j} \mathrm{is} \mathrm{a} \mathrm{nonstream} \mathrm{pixel}$$

(6)

$$\mathrm{or} {\mathrm{B}}_{\mathrm{sum},\mathrm{i}}={ \mathrm{L}}_{\mathrm{sum},\mathrm{i}}{\sum }_{\mathrm{j} \in \left\{\mathrm{cells} \mathrm{to} \mathrm{which} \mathrm{cell} \mathrm{i} \mathrm{pours}\right\}}\mathrm{pij}, \mathrm{if} \mathrm{j} \mathrm{is} \mathrm{a} \mathrm{stream} \mathrm{pixel}$$

(7)

$${\mathrm{B}}_{\mathrm{i}}= \max \left({\mathrm{B}}_{\mathrm{sum},\mathrm{i} }. \frac{{\mathrm{L}}_{\mathrm{i}}}{{\mathrm{L}}_{\mathrm{sum}, \mathrm{i}}}, 0\right)$$

(8)

where L_sum,i is the cumulative upstream recharge, p_ij is the proportion of flow from cell i to j, and the base flow B_i can be directly derived from the proportion of the cumulative base flow leaving cell i, with respect to the recharge available to the upstream cumulative recharge.

2.3.2. Sediment Delivery Ratio Model

InVEST’s SDR model was used to quantify the contribution of sediments to watercourses. The model was derived from studies by Borselli et al. [51] and is based on the quantification of the annual rate of soil loss per pixel and the sediment delivery rate (SDR), which represents the proportion of soil loss that is transported and deposited in watercourses. The model does not consider in-stream processes, assuming that all sediment that reaches the watercourse is transported to the outlet of the hydrographic basin.

In the SDR model, the quantification of the annual soil loss per pixel (ton ha⁻¹·year⁻¹) is estimated using the Revised Universal Soil Loss Equation (RUSLE) (Equation (9)). The determination of the RUSLE parameters is presented in the Supplementary Material.

USLE_i = R_i · K_i · LS_i · C_i · P_i

(9)

where R_i is the rainfall erosivity in pixel i (MJ⋅mm⋅ha⁻¹⋅h⁻¹) (Equation (S2)), K_i is the soil erodibility at pixel i (ton⋅ha⋅h⋅MJ⁻¹⋅ha⁻¹⋅mm⁻¹) [52], LS_i is the slope length-gradient factor at pixel i (dimensional) (Equation (S3)), C_i is the cover-management factor in pixel i (Table S2), and P_i is the support practice factor in pixel i (Table S2).

The sediment delivery rate is calculated as a function of the hydrological connectivity of the basin (Equation (10)). For its determination, it is first necessary to calculate the hydrological connectivity index (IC) (Equation (11)), which describes the hydrological connection between the sources of sediments in the landscape and the watercourses. The higher the pixel’s IC value, the greater is the probability of sediment reaching the watercourse. The IC is based on the relationship between the characteristics of the area upstream of each pixel (D_up), such as land cover, slope, and drainage area (Equation (12)), and the characteristics of the flow path between the pixel and the watercourse (D_dn), such as distance, land cover, and slope (Equation (13)).

$${\mathrm{SDR}}_{\mathrm{i}}=\frac{{\mathrm{SDR}}_{\max }}{1+\exp \left(\frac{{\mathrm{IC}}_0-{\mathrm{IC}}_{\mathrm{i}}}{{\mathrm{k}}_{\mathrm{b}}}\right)}$$

(10)

where SDR_max is the maximum theoretical SDR, set to a mean value of 0.8 [53], and IC₀ and k_b are calibration parameters.

$$\mathrm{IC}={ \log }_{10}\left(\frac{{\mathrm{D}}_{\mathrm{up}}}{{\mathrm{D}}_{\mathrm{dn}}}\right)$$

(11)

$${\mathrm{D}}_{\mathrm{up}}= \overline{\mathrm{C}}\overline{\mathrm{S}}\sqrt{\mathrm{A}}$$

(12)

$${\mathrm{D}}_{\mathrm{dn}}={{\displaystyle \sum }}_{\mathrm{i}}^{}\frac{{\mathrm{d}}_{\mathrm{i}}}{{\mathrm{C}}_{\mathrm{i} }{\mathrm{S}}_{\mathrm{i}}}$$

(13)

where $\overline{\mathrm{C}}$ is the average of the factor C of the upstream catchment area $\overline{\mathrm{S}}$ is the mean of the upstream slope gradient (m. m⁻¹), A is the upstream contribution area (m²), d_i is the length of the flow path along according to the steepest downslope direction (m), C_i is the C factor for each pixel, and S_i is the pixel slope gradient.

Calculation of the exported sediment load per pixel in ton ha⁻¹ year⁻¹ is given by Equation (14). The basin’s total sediment export is calculated by the sum of all pixels constituting the watershed.

$${\mathrm{E}}_{\mathrm{i}}={ \mathrm{USLE}}_{\mathrm{i}}\times { \mathrm{SDR}}_{\mathrm{i}}$$

(14)

2.3.3. Nutrient Delivery Ratio Model

The NDR model was used to quantify annual TN and TP cargo exports. The model maps the nutrient sources along the watershed and their transport to watercourses, enabling quantification of the nutrient retention services provided by vegetation. The model is based on a simple mass balance approach, representing steady-state nutrient flux through empirical relationships. The mass balance is determined based on the nutrient loadings and retention properties of pixels belonging to the same runoff. The algorithm first calculates the annual nutrient load and then the NDR, which corresponds to the proportion of nutrients that reach the watercourse.

Nutrient loads in each pixel are calculated using average TN and TP export coefficients, which vary according to the watershed’s land use. In this study, the coefficients of the Mathematical Model of Correlation between Land Use and Water Quality (MQUAL), version 1.5 [54] were used. The nutrient loads associated with each land-use class were corrected as recommended by Sharp et al. [16] (Equation (S4)). The adjusted values of the nutrient export coefficients are presented in Table S2 in the Supplementary Material.

Then, the model modifies the pixels’ nutrient loads based on the basin’s surface runoff potential (Equation (15)). The model enables analysis of nutrient loads associated with sediments, and that are dissolved, which are then transported through the surface and underground flow, respectively. This study considered only the surface transport of nutrients.

$$\mathrm{modified}.{\mathrm{load}}_{\mathrm{xi}}={ \mathrm{load}}_{\mathrm{xi}}\times {\mathrm{RPI}}_{\mathrm{xi}}$$

(15)

$${\mathrm{RPI}}_{\mathrm{xi}}=\frac{{\mathrm{RP}}_{\mathrm{i}}}{{\mathrm{RP}}_{\mathrm{av}}}$$

(16)

where modified.load_xi is the modified load of nutrients in each pixel, load_xi is the nutrient load on each pixel, RPI_xi is the runoff potential index at pixel i, RP_i is the runoff value of pixel I, and RP_av is the average value of runoff in the watershed. In practice, the RP values are defined through InVEST’s SWY model or through spatialized precipitation data. This study used average annual precipitation data as a proxy to determine the basin’s surface runoff potential.

After quantifying nutrient loads, the model simulates their transport to the watercourse using the NDR factor (Equation (17)), calculated for each pixel based on the IC and retention properties of pixels that belong to the same flow path, a similar approach to the SDR concept [16].

$${\mathrm{NDR}}_{\mathrm{i}}={ \mathrm{NDR}}_{0,\mathrm{i}}{\left(1+\left(\frac{{\mathrm{IC}}_{\mathrm{i}}-{ \mathrm{IC}}_0}{\mathrm{k}}\right) \right)}^{-1}$$

(17)

where NDR_0,i is the nutrient proportion not retained by downstream pixels (see Sharp et [16] for further details), IC_i is the hydrological connectivity index, and IC₀ and k are calibration parameters.

The nutrient load exported per pixel in (kg ha⁻¹ yr⁻¹) is calculated according to Equation (18). The total nutrient export is calculated by the sum of all pixels in the watershed. As InVEST’s nutrient export model only considers nonpoint sources, it was necessary to incorporate TN and TP point sources. Determination of the basin’s point pollution is presented in greater detail in the Supplementary Material.

$${\mathrm{X}}_{\mathrm{expi}}= \mathrm{modified}.{\mathrm{load}}_{\mathrm{xi}}\times {\mathrm{NDR}}_{\mathrm{i}}$$

(18)

2.4. Sensitivity Analysis

Sensitivity analysis was conducted by varying the models’ calibration parameters. For the SWY model, the approach of Hamel et al. [28], in which the α and β parameters were changed one at a time to verify their BF sensitivity, was adopted. With the default value of α = 1/12, the values of β were varied from 0 to 1 in increments of 0.2. The same procedure was performed for α, in which the default value of β = 1 was maintained and the values of α = 1/6 and α = 1/3 were used to analyze their BF response. For the SDR model, parameters k_b and IC₀ (default values: 2 and 0.5, respectively) were modified to assess their influence on annual sediment export. The k_b and IC₀ values were changed from 0.5 to 3 in 0.5 and from 0.25 to 1 in 0.25 increments, respectively. For the NDR model, the parameter k was varied, using the same approach as the SDR model and the parameter critical flow length, that is, the minimum metric distance required for a fragment of a certain type of soil use to retain nutrients maximally. As input data for the NDR model, the critical flow length does not have a default value. The default value of 150 m was considered and the sensitivity of this parameter was set to values of 30, 90, 150, 300, and 500. For all models, the threshold flow accumulation value was set to 1000, the value that best represents the drainage mapped in the basin at a scale of 1:250,000 [55].

2.5. Model Calibration and Validation

This study calibrated and validated the InVEST models to determine the values of the calibration parameters that best represent WES’s spatial and temporal variability. Calibration and temporal validation were performed comparing simulated and observed annual mean values. As the basin’s monitoring stations are sparse, viz., four for streamflow measurement and three for water quality, it was not feasible to use each station’s long-term average values to calculate spatial assessment. Accordingly, to represent spatial variability, the study used all values of the historical series sampled in the basin. WES-specific values were used to normalize the data and eliminate the drainage area effect.

The SWY model was calibrated and validated for the basin’s four sub-basins (SF-I, SF-II, SF-III, and SF-IV), using observed daily and monthly streamflow data [47,48]. The historical series of all sub-basins had faults that varied according to the monitoring station. To fill in the gaps and standardize the 39-year time scale for all stations, the study used a flow regionalization method based on area proportionality, which assumes the existence of a proportional linear relationship between drainage area and flow, a method validated for the region in Anjinho et al. [56,57]. The average annual streamflow (QF + BF) results in mm year⁻¹ simulated by the SWY model were converted to m³ s⁻¹ and compared with observed values. The calibration period was 1981 to 2010 and the validation period was 2011 to 2019.

The SDR and NDR models were calibrated and validated for three sub-basins (WQ-I, WQ-II, and WQ-III), using data from bi-monthly measurements of TS, DS, TN, and TP concentrations. Total suspended solids (TSS) data were obtained indirectly by subtracting TS—DS. Of all the monitoring stations, only WQ-II presents streamflow measurements along with water quality parameters. Thus, for the other sub-basins, the instantaneous streamflow was estimated through regionalization based on the proportionality of the area, using streamflow data from nearby monitoring stations [47]. The annual loads of SST, TN, and TP were estimated based on the concentrations of water quality parameters (mg L⁻¹) and flow (m³ s⁻¹). Historical series of 9 years for sub-basin WQ-I and 19 years for sub-basin WQ-II (12 and 7 years for calibration and validation, respectively) were obtained for both variables. For sub-basin WQ-III, 20 years (12 and 8 years for calibration and validation, respectively) were obtained for TSS, and 30 years (22 and 8 years for calibration and validation, respectively) for TN and TP. Due to the dearth of sampling (merely 6 measurements per year) and the data’s seasonality, the observed loads of TSS, TN, and TP were processed to represent average annual values. For TN and TP, the outliers were removed using the relation lower outlier = Q1 − (1.5 × IQR) and higher outlier = Q3 + (1.5 × IQR), where Q1 is quartile 1, Q3 is quartile 3, and IQR is interquartile variation (IQR = Q3 − Q1). For TSS, in which a greater effect of seasonality was observed, only the values that were between Q1 to Q3 were considered.

2.6. Model Performance Analysis

The performance of the models was evaluated by comparison with data observed in the monitoring stations, considering annual average values. The analysis was conducted to evaluate WES temporal and spatial performance. Three statistical indicators were used: the Nash and Sutcliffe coefficient [58], percentage bias (PBIAS), and the coefficient of determination (R²). The classification suggested by Rauf et al. [59] was used to interpret the results (

Table 1. Classification of statistical performance indicators.

Classes	NSE	R2	PBIAS
Very good	0.75 > NSE ≤ 1.00	0.85 > R2 ≤ 1.00	<|5|
Good	0.65 > NSE ≤ 0.75	0.70 > R2 ≤ 0.85	|5| − |10|
Satisfactory	0.50 > NSE ≤ 0.65	0.60 > R2 ≤ 0.70	|10| − |15|
Acceptable	0.40 > NSE ≤ 0.50	0.40 > R2 ≤ 0.60
Unsatisfactory	NSE ≤ 0.40	≤0.40	≥|15|

).

3. Results

3.1. WES Models Sensitivity Analysis

The sensitivity analysis of the SWY model parameters indicated that β is the most sensitive and has the greatest influence on BF variation (Figure 3). Increased β values decreased BF in all sub-basins, primarily, in the range of values that varies between 0 ≤ β < 0.4. Compared to the default values, reduction in the value of β to 0.2 increases BF 26% on average, with the greatest variation (39%) in sub-basin SF-IV. For α, however, the trends varied. BF decreased by an average of 3% to α = 1/6 and increased by an average of 5% to α = 1/3, chiefly in sub-basin SF-IV, where the average increase was approximately 10%.

Sediment export in the SDR model showed greater sensitivity for parameter k_b, whose increase is associated with greater sediment export to all sub-basins (Figure 3). For k_b = 0.5, the exported sediment value is on average 100% lower than the default value, and for k_b = 3, the value is on average 107% higher. In absolute values (t yr⁻¹), the highest sensitivity occurs between 0.5 ≤ k_b ≤ 1, where the sediment exported for k = 1 is 34 times greater. Values between 1.5 ≤ k_b ≤ 2.5 tend to vary sediment export by approximately −50% and +50%, respectively. Changes in IC0 values showed lower sensitivity and opposite behavior compared to k_b variations. The increase in IC₀ values is related to a lower sediment export. On average, an increase of 0.25 in IC₀ values tends to decrease the value of exported sediment by 11%. Compared to the default value, the value of IC₀ = 1 provided an increase of 20% of the sediment exported to all sub-basins, while the value of IC₀ = 0.25 decreased by 12%.

The sensitivity of the parameters k and critical flow length showed the same behavior for TP and TN (Figure 3). In general, nutrient exports were more sensitive to changes in critical flow length values. This study considered the value of 150 m as the default for this parameter, decreasing this value to 30 m provided a reduction of approximately 12% in the value of exported nutrients, while the value of 500 m increased by almost 25%. The parameter k was more sensitive between 0.5 ≤ k ≤ 1 and its effect is more expressive in the decrease of nutrient export than in the increase. Relative to the default value (k = 2), decreasing to k = 0.5 decreased nutrient export by 15%, while changing to k = 3 showed a slight increase of almost 3%.

3.2. Calibration and Validation of WES Temporal Variability

The SWY, SDR, and NDR models perform differently in the spatial and temporal scales. Figure 4, Figure 5 and Figure 6 show the simulated and observed WES values over time and the models’ adjusted parameter values.

Table 2. InVEST WES model performance in light of temporal variability (* inversely proportional relationships, negative slope).

			Calibration			Validation
Monitoring Station	Model	Indicator	R2	PBIAS (%)	NSE	R2	PBIAS (%)	NSE
SF-I	SWY	Streamflow	0.4	−0.6	−0.6	0.5	8.1	0.6
SF-II	SWY	Streamflow	0.5	6.6	−0.6	0.8	30.3	0.6
SF-III	SWY	Streamflow	0.5	8.8	−4.6	0.7	0.9	0.5
SF-IV	SWY	Streamflow	0.6	−3.4	−2	0.7	−8.5	0.6
WQ-I	SDR	Exported sediment	0.3	0.2	0.2	-	-	-
WQ-II	SDR	Exported sediment	0	0.6	−0.1	0.2	−14.6	−0.1
WQ-III	SDR	Exported sediment	0	6.6	−0.3	0.1 *	5.6	−0.4
WQ-I	NDR	Exported TP	0.2	50	−4.5	-	-	-
WQ-II	NDR	Exported TP	0.6 *	−23	−1.1	0.4	1.7	0.2
WQ-III	NDR	Exported TP	0.3	−10	0.1	0.3 *	2.2	−0.3
WQ-I	NDR	Exported TN	0.3	−18.6	−0.8	−	−	−
WQ-II	NDR	Exported TN	0.4 *	2	−0.2	0.4 *	12.3	−0.6
WQ-III	NDR	Exported TN	0.1	27	−1.5	0.6 *	44	−3.4

presents performance indicators. The best calibration and validation performance was observed for the SWY model. A visual comparison of the simulated and observed annual streamflow indicates a good fit for all sub-basins (Figure 4). Although the best performance involved values close to the long-term annual averages, the model was able to represent streamflow variation over time. R² indicated an acceptable performance for the calibration period, but NSE was unsatisfactory for all sub-basins. PBIAS indicated that the SWY model underestimated annual streamflow in sub-basins SF-I and SF-IV and overestimated it in sub-basins SF-II and SF-III, with very good and good performance, respectively. For the validation period, R² values indicated satisfactory to acceptable performances. NSE values indicated satisfactory performance and PBIAS values presented different values for each sub-basin, with unsatisfactory performance observed in sub-basin SF-II and very good performance in sub-basin SF-III.

The SDR model effectively represented years with values close to the long-term average, but did not perform well in representing the annual variability of exported sediment (Figure 5). The results indicated large discrepancies for specific years, but those between the observed and simulated values were not systemic. NSE and R² values for the calibration and validation period reflected the SDR model’s limited capacity to represent the annual variability of exported sediment, indicating inconsistent performance (Table 2). PBIAS values showed that the SDR model overestimated the sediment load exported annually in the calibration period, with very good performance in sub-basins WQ-I and WQ-II and good performance for sub-basin WQ-III. PBIAS found varying trends in the sub-basins for the validation period, underestimating the sediment load exported in sub-basin WQ-II and overestimating it in sub-basin WQ-III, with satisfactory and good performance, respectively (Table 2).

Nutrients simulated by the NDR model present higher performance for values close to the long-term mean of the observed data. However, unlike the SWY and SDR models, the NDR model is unable to represent the annual variability of nutrient loads (Figure 6). For TP and TN, values tend to remain close to the average, showing a slight variation over time. Concerning the NDR model’s performance, the adjustment of the calibration parameters showed high variability among sub-basins. With regard to the values of the annual TP loads for the calibration and validation period, the R² and NSE indicators indicated unsatisfactory performance for all sub-basins (Table 2). Although R² was increased for sub-basin WQ-II (R² = 0.6) in the calibration period, the analysis found an inversely proportional relationship between the simulated and observed data, which indicates the model’s inconsistent performance. PBIAS, however, indicated satisfactory performance in the calibration period for sub-basin WQ-III and very good performance for the validation period in sub-basins WQ-II and WQ-III (Table 2). The performance indicators obtained for modeled TN annual loads showed the same trend as for modeled TP annual loads. The PBIAS indicator indicated very good and satisfactory performances only in sub-basin WQ-II for the calibration and validation period, respectively (Table 2).

3.3. WES Spatial Performance Analysis

SWY model results demonstrated a reasonable fit with the observed data (

Table 3. InVEST WES model performance considering spatial variability.

Model	R2	PBIAS	NSE
SWY	0.47	3.50	−0.39
SDR	−0.10	−6.25	−0.13
NDR-TN	0.30	11.00	0.13
NDR-TP	0.00	−3.31	−0.23

). R² indicated acceptable performance, and the PBIAS and NSE indicators demonstrated a very good and unsatisfactory performance, respectively. For the SDR, NDR-TN, and NDR-TP models, PBIAS indicated a performance ranging from satisfactory to very good, while the other indicators showed unsatisfactory performance.

The study’s results indicate that, in general, the models were able to simulate the long-term mean values observed in the sub-basins (

Table 4. Long-term values of simulated and observed data of variables QF, BF, streamflow, exported sediment, TN, and TP.

Sub-Basin	Period	QF (m³ s−1)	BF (m³ s−1)	Simulated Streaflow (m³ s−1)	Observed Streamflow (m³ s−1)	Discrepancy (%)
SF-I	1981–2019	1.2	2.6	3.8	3.8	1
SF-II	1981–2019	3.1	9.4	12.5	11.3	11
SF-III	1981–2019	8.9	18.0	26.9	25.1	7
SF-IV	1981–2019	13.1	32.6	45.8	47.9	−5
Sub-basin	Period	Simulated exported sediment (kt.yr−1)	Observed exported sediment (kt.yr−1)			Discrepancy(%)
WQ-I	2011–2019	8.350	7.759			8
WQ-II	2001–2019	14.843	15.533			−4
WQ-IIII	2000–2019	27.557	31.757			−13
Sub-basin	Period	Nonpoint TN load (t.yr1)	Point TN load (t.yr1)	Simulated TN (t.yr1)	Observed TN (t.yr1)	Discrepancy(%)
WQ-I	2011–2019	136.4	280.4	416.8	512.0	−19
WQ-II	2001–2019	242.3	331.4	573.6	542.4	6
WQ-IIII	2000–2019	485.9	707.2	1193.0	900.0	32
Sub-basin	Period	Nonpoint TP load (t.yr1)	Point TP load (t.yr1)	Simulated TP (t.yr1)	Observed TP (t.yr1)	Discrepancy(%)
WQ-I	2011–2019	16.2	27.0	43.2	28.8	50
WQ-II	2001–2019	29.6	32.6	62.2	73.2	−15
WQ-IIII	2000–2019	63.0	71.7	134.8	143.9	−6

). The largest discrepancy occurred in sub-basin WQ-I for the NDR-TP model, and the smallest discrepancy was observed in sub-basin SF-I for the SWY model.

4. Discussion

4.1. WES Sensitivity Analysis

Unlike other studies that analyzed the sensitivity of calibration parameters and input data [27,60,61], this study examined just the calibration parameters’ sensitivity, choosing not to analyze the sensitivity of the input data because the vast bulk came from monitoring stations, which tend to have higher quality than other databases. When feasible, empirical values were also used for some input parameters, such as soil erodibility and crop and nutrient export coefficients, which best represent the characteristics of the studied area.

With regard to the SWY model, the study indicated the β parameter as the most sensitive, having a significant effect on flow routing and, consequently, on BF values. Hamel et al. [28] also found greater sensitivity for β, with BF varying by nearly 50% in some sub-basins. For α = 1/12 and α = 1/6, the decreasing BF was similar to Hamel et al. [28], but for alpha = 1/3, there was a BF increase (Figure 3), associated with the determination of AET_i,m, which is conditioned by PET or water availability, according to Equation (4). The parameters α and β are associated with the fraction of the annual recharge of pixels upstream (upslope) available for evapotranspiration of a pixel downstream (downslope) [16]. The configuration of alpha = 1/3 and beta = 1 (default value) renders the potential available water (L_sum.avail,i), very high, requiring the algorithm to use PET_i,m to calculate AET_i,m (Equation (4)).

The sensitivity of the SDR model parameters k and IC₀ followed that of Hamel et al. [25], as anticipated in light of the model’s configuration [16]. Although both parameters characterize the relationship between the connectivity index and sediment delivery rate, the sensitivity of k was greater than that of IC₀, with higher values associated with a higher sediment delivery rate, as predicted by Vigiak et al. [53], suggesting that k variations are more effective for data calibration and validation. The result of the sensitivity analysis of the IC₀ parameter appears more sensitive than that in Hamel et al. [25], although it showed the same linear decrease in exported sediment.

For the NDR model, the sensitivity of the parameter k occurred as anticipated, given the relationship between nutrient export and k is generally exponential when k is less than nine and linear when it is greater [62]. As noted in other studies [27,62], the results of this research showed that the parameter k is most sensitive between the values of 0 ≤ k ≤ 1, in which the projected curve presents sigmoidal behavior, that is, an expressive increase at the onset and a tendency to stabilize as k values increase. In this study, the effect of k in all sub-basins was similar (Figure 3) as a result of basin’s hydrogeological similarity. However, as demonstrated by Redhead et al. [27], the response in k variations is specific to each location, and thus may differ from the results herein. Variations in critical flow length values indicate a directly proportional linear relationship with the exported nutrient load. This was also expected, since it is assumed that the greater the distance for a land-use fragment to reach its maximum retention capacity, the greater the probability of the nutrient load reaching the watershed [16].

4.2. InVEST WES Model Performance

The study’s sensitivity analysis increases our understanding of the effect of calibration parameters on the modeled results and how to optimize their representation of the observed data, which its performance analysis indicates the InVEST WES models are reasonably effective in simulating, albeit with high variability according to the model and analysis scale. Using the SWAT model and the Generalized Watershed Loading Function, Santos et al. [63] modeled the streamflow and sediment and nutrient loads in the upper part of the basin with results analogous to those herein, with both indicating improved performance in simulating streamflow.

With regard to the InVEST WES models’ calibration and temporal validation, the best performance was observed for the SWY model. Despite some indicators of unsatisfactory performance, particularly in the calibration period for NSE, the results show that the model was able to represent the temporal variability of the average annual streamflow in the sub-basins for the calibration and validation period (Figure 4 and Table 2). The effective performance of the SWY model may be associated with its high sensitivity to variations in precipitation, as noted in the literature [29,61]. Likewise, the inferior performance of the NDR and SDR models may be related to their lower sensitivity [27,64,65]. Insensitivity to precipitation in UK watersheds [27] is inherent in the NDR model’s configuration (Equation (16)), where its sole effect is to modify the nutrient load’s impact on surface runoff, based on an index that relates pixel precipitation to the basin’s average rainfall. That is, the effect of precipitation is more closely associated with its configuration than with its intensity. In the SDR model, as opposed to the NDR, precipitation has a greater expressive influence on sediment export (Figure 5) because rainfall erosivity is key to triggering erosion (Equation (9)). Some studies have indicated high sensitivity to erosivity in the SDR model [25,66]. Although this study used annual rainfall and land-use data as model inputs, temporal variability was not well captured by the SDR and NDR models.

The spatial performance of InVEST’s WES models is often evaluated using long-term average values. Due to the dearth of monitoring stations, this study included all values in the historical series sampled in the basin to calculate the statistical indicators and evaluate the spatial performance of the models. Although not detailed herein, it should be noted that an initial analysis was conducted, using WES values not normalized by the drainage area, which deemed the performance of all models good. However, as the study’s objective was to examine the models’ capacity to forecast WES behavior in the basin, taking into account the variability of land use patterns, specific WES values were used, and the models’ spatial performances were significantly reduced (Table 3).

When compared with observed long-term mean values, the results were similar to those in the literature [13,24,25,27,28], indicating that the InVEST WES models yield good results in terms of relative magnitude. The SWY model performed best, with simulated long-term average flow deviations of less than 10% from observed values. These results confirm the findings of other studies that have shown reasonable performance for the SWY model [13,28,29]. The long-term average loads of exported sediments simulated by the SDR model also showed limited discrepancies with observed values, corroborating studies that have featured its potential to simulate sediment dynamics [24,25,66]. The greatest deviations in NDR model simulations were found in sub-basins WQ-I and WQ-III, which overestimated TN and TP exports by 50% and 32%, respectively. Redhead et al. [27] obtained similar results in UK watersheds, with the vast majority showing deviations from observed values for TN and TP of 65% and 44%, respectively.

In calibrating the models, it could be seen that that the values assigned to the calibration parameters generally reflected the hydrogeological characteristics of the basin, which is located in a flat region, in which some 80% of the relief is composed of smoothly undulating hills with low drainage density [43], with a slope around 8%. The basin’s soils are primarily sandy and deep (~75%), as is characteristic of Brazil’s Cerrado biome. Thus, the basin has a low potential to generate surface runoff and, consequently, to transport nutrients and of sediments to watercourses. The values assigned to the calibration parameters reflected these characteristics.

In the SWY model, the values β = 0.2 and α = 1/12 presented the best relationship with observed values and increased the BF share of water yield, conditions that occur in regions with flat, smooth topography. Their determination is consistent for the basin as it is an important recharge area for the Guarani Aquifer, one of the largest underground freshwater reservoirs in the world [41,42]. The values assigned to the SDR model parameters (k = 0.8 and IC0 = 0.45) also reflect the basin’s environmental characteristics, characterizing environments with low potential for sediment transport, although their soils are highly susceptible to erosion [57,67]. The NDR model-TP was the only one in which the default value of the calibration parameter was maintained (k = 2), although the transport of nutrients was reduced through the critical flow length parameter.

4.3. InVEST WES Model Potentials and Limitations

InVEST’s WES models present simple approaches and demand few input data compared to other physical-based hydrological models, such as SWAT. These characteristics facilitate decision making in non-instrumented basins or regions where data is scarce, as demonstrated by Benra et al. [13] and corroborated by this study. InVEST’s WES models enable spatial understanding since they are distributed at pixel scale, and the WES value in each pixel varies depending on its location, considering such key factors as land use, vegetation, relief, and soil types, which remain a major challenge for ecohydrological models [68,69]. Such characteristics make these models important tools for environmental planning, enabling an analysis of the effects of diverse future scenarios, such as climate and land-use changes, the primary drivers of WES degradation [4,5]. While WES model performance differs for each model and temporal and spatial scale, in general, the study results indicate the models’ capacity to represent long-term average values, and, in the case of the SWY model, to simulate annual average values, enabling a representation of the variability of extreme annual values (see 1983 in Figure 4).

Like any other model, InVEST’s WES models have advantages and limitations. Accordingly, further research to enhance them and increase their transparency is in order. As noted in the InVEST manual, the models use simple equations and limited parameters to estimate water ecosystems services [16]. An analysis of the results of this study confirm the limitations noted in the InVEST user manual and in other studies [16,25,27]. The principal limitation of the SWY model is its inability to accurately represent water flow, that is, surface runoff and base flow. The discrepancies arise from its simulation of surface runoff, based on an adaptation of the SCS curve number method, which disregards topographical effects, and of base flow, whose simplified routing equations can generate inconsistent results in terms of absolute numbers [70]. It should also be noted the model fails to simulate the annual variability of water yield in the watershed. Although it enables evaluation of the monthly surface runoff, this limits its use to contexts that demand data on a more refined time scale, as in the simulation of reference streamflow to calculate water availability [71].

The effect of simplification is even more significant in the SDR and NDR models, which are even more sensitive to input data. Uncertainties related to input parameters, such as R and K from RUSLE, and export and nutrient retention coefficients can have considerable impact on the models’ simulations [16,25]. In many regions, there is no information available on these empirically determined parameters, necessitating the use of data from other locations, which may have quite different climatic and hydrogeological characteristics.

Other limitations associated with the generalization of complex physical processes relate to nutrient and sediment dynamics. In the SDR model, sediment yield is solely determined by RUSLE, disregarding such sources as the gully, stream bank, and mass erosion, which are useful in comparing simulations with observations [25]. The determination of TN and TP loads in the NDR model disregards the cycle of each nutrient, quantifying both loads in the same way through the relationship between land-use and export coefficients, a considerable abstraction given that their environmental dynamics are quite different [72]. Neither the SDR nor NDR model take into account the dynamics of nutrients and sediments in watercourses, assuming that the entire load that enters the watercourses is exported to the watershed. The NDR model does not take into account point loads, which can adversely impact urbanized basins, especially in areas where basic sanitation is deficient. In the basin, for example, where urban areas represent only 4% of the area, points contribute, on average, up to 60% of the total loads.

4.4. Potential, Uncertainties, and Limitations

Unlike most studies that analyzed the spatial performance of InVEST’s WES models, considering long-term average values [13,24,25,26,27,28], this study analyzes their performance, considering WES spatial and temporal variability. Of the studies reviewed, only Lu et al. [29] analyzed the temporal performance of the SWY model using Pearson’s correlation, albeit using a brief historical series of observed flow data. This study uses historic series of land-use and meteorological data measured at monitoring stations to simulate the temporal variability of WES, considering climate and land use, which exert significant influence on WES dynamics in hydrographic basins.

The principal uncertainties pertaining to the study’s calibration method pertain to the use of historical series of observed data. The daily and monthly streamflow data, for example, had numerous flaws, and several stations lacked data for the complete series (1981–2019). With regard to TN, TP, and TSS water quality data, it should be noted that CETESB samples are only taken every two months. In addition, many CETESB monitoring stations do not perform streamflow measurements at the points where water quality samples are taken, which makes determining nutrient and sediment loads needed for comparison with model results problematic. An attempt was made to circumvent flow data deficiencies using hydrological regionalization, and the effects of limited annual sampling and, consequently, the seasonality of the water quality data, were reduced by removing outliers (see Section 2.5), since bi-monthly measurements do not represent the annual TN, TP, and TSS averages. Thus, despite strenuous efforts to calibrate and validate the WES models, uncertainties and limitations remain.

In addition to such challenges with observed data, there are uncertainties related to the empirical parameters used in the models. The study sought to use regional data for all parameters, since they best represent the environmental characteristics of the basin studied (see Supplementary Material). However, for some parameters, such as CN and nutrient retention efficiency coefficients, for example, it was not feasible to use empirical values, which can contribute to uncertainty in model results. For the NDR model, uncertainties include quantification of TN and TP point loads, estimated from on the per capita load of domestic sewage, using export coefficients, population data, and nutrient removal coefficients, and the values are uncertain due to the scarcity of data on basic sanitation in the basin’s municipalities.

5. Conclusions

This study investigated the performance of InVEST 3.9.0 WES models, analyzing their ability to represent the spatial and temporal variability of observed values. The results show a high variability in their performances, according to the type of water ecosystem services simulated and scales of analysis. The best performance was observed for the SWY model, which effectively represented the spatial and temporal variability of the average annual streamflow in the analyzed sub-basins. The SDR and NDR models were unable to represent the temporal variability of the exported loads of sediments and nutrients. The inferior performance of these models may be associated with lower sensitivity to precipitation variability. The NDR model, in particular, appears insensitive to annual variations in precipitation, which impact the annual load of nutrients exported in the basin.

The spatial performance considering specific WES values was low for most statistical indicators. In general, with the exception of some sub-basins, all models showed good performance in simulating long-term mean values, which corroborates the results of other studies.

Properly calibrating models’ parameters is essential to enhancing their performance, and the values the study assigned to them reflected the basin’s hydrogeological characteristics. Uncertainties regarding the methodology used for calibration and validation relate to the quality of the historical series of observed data of daily and monthly streamflow and inadequate water quality sampling data, as well as the failure to measure streamflow data. While the study attempted to circumvent such obstacles, uncertainties remain, and the results should be interpreted with circumspection. Furthermore, some findings may be solely empirical, as the study area lacks broad hydrogeological variability.

Despite the uncertainties and empirical nature of this study, its results and discussions should aid researchers and other users of InVEST models to inform and make decisions and formulate policies. To minimize uncertainties and enhance results, the use of the complete historical series of streamflow and, particularly of water quality parameters is recommended to better reflect the basin’s hydrological characteristics and correspond to simulated data. In as much as the performance of the models may vary according to the hydrogeological nature of the region in which they are being applied, further research in varied basins in diverse regions should be conducted to complement the findings detailed herein.