WRF-Hydro for Streamflow Simulation in the MATOPIBA Region within the Tocantins/Araguaia River Basin—Brazil: Implications for Water Resource Management

Abstract

The effective management of water resources in regions with a high potential for water resources, such as the Tocantins/Araguaia Basin in Brazil, is crucial in the face of current climate change and urban and agricultural expansion. In this context, this study evaluates the WRF-Hydro hydrological model to simulate the flow of the Manuel Alves Pequeno, Vermelho, and Manuel Alves Grande rivers in the MATOPIBA region (encompassing areas from the states of Maranhão, Tocantins, Piauí, and Bahia), an agricultural frontier and the most key area in terms of grain production in Brazil. The aim is to analyze the hydrological parameters of soil infiltration, surface retention depth, land surface roughness, and Manning’s channel roughness. The simulations are conducted at a spatial resolution of 3 km with a channel network of 100 m, covering a period of heavy rainfall from 13 March to 1 June 2018. For model validation, observational data from three river gauge stations of the National Water and Sanitation Agency are used, with assessments based on the Nash-Sutcliffe efficiency index, standard deviation of observations, root mean square error, percentage bias, and correlation coefficient, resulting in values of 0.69, 0.56, 4.99, and 0.83, respectively. In particular, the adjustment of the infiltration factor and surface roughness parameter has a greater contribution to improving the statistical results than the adjustment of the other two hydrological parameters. Additionally, the quality of discharge simulation at each river gauge station is correlated with the temporal distribution of simulated precipitation compared to observed data in the drainage network. Highlighting WRF-Hydro’s potential as a fine-scale model easily coupled with numerical weather prediction, this study significantly advances regional river dynamics evaluation, crucial for strategic water resource management.

1. Introduction

Ensuring water supply is an escalating concern in the current global context, mainly due to the exacerbation of water scarcity in specific regions. In Brazil, the water crisis, primarily faced in major urban centers, arises from urban expansion driven by population growth, increasingly intermittent precipitation patterns due to climate change, and inefficient water resource management. Within this framework, the rising demand for water poses a significant challenge for both human consumption and activities related to agriculture, livestock, industries, and the energy sector [1].

Brazil is among the countries with the highest availability of freshwater resources globally [2] and ranks as the fourth-largest grain producer globally, surpassed only by China, the United States, and India. It leads in producing and exporting soybeans, accounting for 50% of global trade [3]. A region significantly contributing to this productivity is the new national agricultural frontier known as MATOPIBA, depicted in Figure 1, an acronym representing the initials of the states it encompasses: Maranhão, Tocantins, Piauí, and Bahia. This area has gained national prominence due to its robust potential for grain and fiber production [4]. Within this region, three river basins are present: the São Francisco River basin, the Atlantic basin—North/Northeast section, and the Tocantins/Araguaia basin, which is the focus of this article, specifically the Tocantins sub-basin between the Sono and Araguaia rivers. The Tocantins-Araguaia river basin, the largest river basin entirely within the national territory, spanning the northern and central-western regions of the country, covers an area of approximately 918.822 km² and occupies 11% of the national territory, including the states of Goiás, Tocantins, Pará, Maranhão, Mato Grosso, and the Federal District [5]. The region has abundant water resources supporting various activities, ranging from irrigation to energy production, mining, and industry, emphasizing processing agricultural products such as soybeans, rice, corn, and beef. Municipalities near the Araguaia River exhibit the highest water demands for irrigation [6]. Even though Brazil has a large water source, the growing water demand, both to satisfy human needs and for agriculture and industry, is influenced by temporal and spatial variability, imposing limits on the availability of water in some areas of the country, which presents a major challenge for water resource planning [7]. The complexity of this scenario underscores the importance of studies aimed at investigating the space-time variability of hydrometeorological variables in watersheds to address and reverse this situation. In this context, hydrological and climatic models are tools that can assist decision-makers in water resource management, promoting water security for present and future generations.

Numerous studies focus on addressing issues in the areas of resource conservation and agricultural productivity, showing interest in the complementary modeling of hydrological and cropping systems [8]. The potential of the LASH (Lavras Simulation of the Hydrology) model was assessed by [9] to simulate hydrological impacts in response to land-use change scenarios, following the most significant trends in agribusiness in an Amazon Basin in Brazil. In contrast, [10] evaluated the effects of agricultural expansion on water availability in a region of western Bahia, located in Brazil and part of the MATOPIBA region, through hydrological modeling by assessing river flows using the coupled INLAND (Integrated Model of Land Surface Processes) and THMB (Terrestrial Hydrology Model with Biogeochemistry).

Hydrometeorological studies, particularly in unexplored regions like MATOPIBA in Brazil, are poised to benefit substantially from sophisticated models like the Weather Research and Forecasting Model (WRF-Hydro). WRF-Hydro, a relatively recent entrant since 2013, has already carved a niche for itself due to its robust, integrated architecture and adaptive capabilities across varied geographic and climatic landscapes, which have been corroborated by extensive validations across different global regions [11].

This research leverages the WRF-Hydro model’s comprehensive approach, seamlessly integrating atmospheric and terrestrial components for a nuanced simulation of water cycles, energy transitions, and hydrometeorological variables. Such a comprehensive perspective is crucial for effectively addressing ambiguities related to the spatial allocation of precipitation intensities and enabling the harmonization of high-definition atmospheric models with high-precision hydrological models [10,12,13,14,15].

A critical aspect of our study involves delving into four pivotal parameters: the infiltration factor (REFKDT), surface retention depth (RETDEPRT), surface roughness (OVROUGHRT), and channel roughness (MannN). Prior research has underscored the significance of these parameters in the nuanced calibration and optimization of the model to capture hydrological responses accurately, especially in scenarios characterized by extreme rainfall and flooding events [16,17,18,19,20]. While previous studies have conducted numerous simulations in various global regions, the utilization of the WRF-Hydro model in the MATOPIBA region of Brazil remains largely unexplored. In the Tocantins-Araguaia River Basin, accurate modeling of natural streamflow provides crucial scientific foundations for improving reservoir management and supporting agricultural planning. Despite the significance of these endeavors, the application of such simulations in this specific region has been notably absent. Consequently, this study seeks to employ the WRF-Hydro model to simulate river streamflow, explore the model’s sensitivity to critical parameters influencing its behavior, and evaluate the model’s performance using calibrated parameters. The findings of this research could serve as valuable reference points for further investigations within the study area and neighboring regions.

2. Materials and Methods

2.1. Study Area and Data

The study area is located in one of Brazil’s major agricultural production hubs, known as MATOPIBA, a region comprising 337 municipalities across the states of Maranhão (MA), Tocantins (TO), Piauí (PI), and Bahia (BA). It spans an area of 731.700 km² and is intersected by three river basins: the Tocantins/Araguaia basin (occupying 43% of the area), the Atlantic basin—North/Northeast section (with 40%), and the São Francisco River basin (with 17%) [21].

The climate in MATOPIBA is humid tropical with a dry winter, with average monthly temperatures ranging from 25 to 27 °C throughout the year and an annual average precipitation between 800 and 2000 mm, distributed across two well-defined seasons: the dry season from May to September and the wet season from October to April. The dominant vegetation in the region is the Cerrado biome, which covers 91% of the area. Additionally, the region encompasses remnants of the Amazon (7%) and areas of the Caatinga (2%) [22].

This hydrological study is conducted in the Tocantins-Araguaia river basin, specifically in one of its seven sub-basins belonging to MATOPIBA, the Tocantins sub-basin between the Sono and Araguaia rivers. This area significantly influences the regional drainage of the Tocantins-Araguaia and São Francisco river basins [21,23,24]. The Araguaia River is the main tributary of the Tocantins River, and the Sono River is another significant tributary in the region. The confluence points of the Araguaia and Sono rivers with the Tocantins River essentially demarcate the northern and southern boundaries of the sub-basin.

The rivers analyzed by WRF-Hydro are the Manuel Alves Pequeno, Vermelho, and Manuel Alves Grande rivers, for which flow data are measured at the Itacajá, Jacaré, and Goiatins stations, respectively. These rivers are tributaries of the Tocantins River, all of them with drainage networks entirely within the Tocantins sub-basin, between the Sono and Araguaia rivers. Figure 1 shows the location of the MATOPIBA region and the study sub-basin with its drainage network.

The simulated discharge data is assessed through a comparison with the actual daily discharge data recorded daily by the river gauging stations operated by the National Water and Sanitation Agency (ANA), and these data are available on their official website (www.snirh.gov.br (accessed on 16 September 2023)) [25].

Although the primary aim of this study is the analysis of river flow, the presentation of simulated and observed precipitation values is also included. This is because the hydrological module of WRF-Hydro in this work is fully coupled with the WRF atmospheric model, indicating that it simulates river discharge data based on the simulated data from the meteorological module, including precipitation, which directly affects the discharge results. Therefore, this inclusion allows us to check whether any potential disparities between the simulated and observed flows are a result of likely differences between simulated and observed precipitation.

The observed precipitation data are obtained from the daily records provided by the Center for Weather Forecasting and Climate Studies (CPTEC) at the National Institute for Space Research (INPE) through the MERGE product. This product combines observed precipitation from meteorological stations with satellite precipitation estimates. The MERGE product data can be accessed on the CPTEC website (http://ftp.cptec.inpe.br/modelos/tempo/MERGE/GPM/DAILY/ (accessed on 10 July 2022)) [26,27,28]. This technique provides precipitation values at regularly spaced grid points, similar to the WRF, enabling an analysis that does not rely solely on individual data points but considers any part of the grid within the study domain.

The simulated and observed precipitation data to be presented are the arithmetic mean of the values at the WRF and MERGE grid points in the region influencing the channel flow under study—a region spanning from the river’s headwaters to its mouth. Figure 2 displays the MERGE grid points in the drainage network of the Manuel Alves Grande River, that is, in the area where surface runoff influences the river’s discharge, as well as the location of the hydrological stations Itacajá, Jacaré, and Goiatins, with altitudes of 250 m, 200 m, and 185 m, respectively.

The main issue when forcing such hydrological models with atmospheric model data remains the limited accuracy of simulated precipitation. Therefore, the temporal accuracy of precipitation is analyzed by calculating the relative error (ER) between the simulated and observed precipitation, as defined by Equation (1) [17]:

$$ER=\frac{P-Q}{Q}\times 100$$

(1)

where P is the simulated value, calculated as the average precipitation of all grid points within the drainage network study area of the respective river, and Q is the observed value, also calculated as the arithmetic mean of the precipitation values presented in the MERGE product in the same region as the simulated values.

2.2. Coupled Atmospheric and Hydrological Modeling System

2.2.1. Atmospheric Model (WRF)

The WRF-Hydro was employed in its fully coupled mode, wherein the atmospheric and hydrological models run simultaneously. The atmospheric model is WRF version 3.9.1.1 and was configured with two nested grids, as shown in Figure 1. The outer domain (D01) covers the central part of the MATOPIBA region, has a resolution of 9 km, and consists of 60 grid points in the west-east direction and 66 in the south-north direction. The inner domain (D02), where the hydrological analysis was conducted and encompasses the entire drainage network of the observed rivers, is located in the Tocantins sub-basin between the Sono and Araguaia rivers. It has a spatial resolution of 3 km and comprises 64 grid points in the west-east direction and 67 in the south-north direction. Fifty vertical levels were used, with the model’s top defined at 50 hPa. Figure 1 displays the location and distribution of the domains in the model and a zoomed-in view of D02, including the drainage network.

The meteorological simulation process was carried out using data derived from the United States scientific agency NCEP’s (National Centers for Environmental Prediction) global GDAS-FNL (Global Data Assimilation System—Final Analysis) system, with a horizontal resolution of 0.25° × 0.25°, 32 vertical levels, and a temporal resolution of 6 h [29]. Land use and cover data were provided by MapBiomas, a multi-institutional initiative aimed at generating annual land use and cover data for all of Brazil through an automated classification process applied to satellite data. A historical series from 1985 to 2019 of land use covering 27 categories was developed based on the LANDSAT collection (with a spatial resolution of 30 m) and is available on their official website (https://mapbiomas.org/ (accessed on 15 July 2022)) [30].

The simulations were conducted over a period encompassing occurrences of intense rainfall from 1 December 2017 to 1 June 2018, such that the first 102 days were considered the spin-up period. The physical parameterizations were selected based on precipitation simulation studies using WRF and WRF-Hydro in regions with a climate similar to that of the Tocantins-Araguaia river basin, leading to the set of parameterizations outlined in

Table 1. Main physical parameterizations used in the WRF model in this study.

Parameters	Scheme	Reference
Microphysics scheme	WSM6	[35]
Cumulus convection	Grell-3D	[36]
Longwave radiation and Shortwave radiation	RRTMG	[37]
Planetary boundary layer	MYNN 2.5	[38]
Surface layer	MYNN	[38]
Land surface scheme	Noah LSM	[39]

[31,32,33,34,35,36,37,38,39].

2.2.2. Hydrological Model (WRF-Hydro)

The WRF-Hydro version is 5.2.1, fully coupled with WRF in the D02 domain. Geoprocessing tools were used in this domain to resize the grid from 3 km to 100 m, aiming to create input data for WRF-Hydro related to surface water, groundwater, and channel flows. These data were obtained with the assistance of the WRF Hydro GIS Pre-Processing Toolkit v. 5.1, developed by NCAR, for use in the GIS environment [40]. This toolkit utilized the WRF Preprocessing System input files and the Digital Elevation Model (DEM) to generate high-resolution fields on routing grids, including flow direction, subsurface flow, and channel routing processes. These fields were used as input data in the WRF-Hydro model. A summary of the data used in the model, along with their interactions and the observational data used for validation, is displayed in Figure 3.

The primary hydrological settings include activating the surface and subsurface runoff modules for all domains, while the channel routing module is activated only in the study sub-basin (D02). The bucket model [41,42,43] is employed for calculating the base flow of the basin. The Noah Land Surface Model (LSM) is the terrestrial surface model in the coupled WRF/WRF-Hydro simulations. This model is responsible for the hydrological processes of the column (i.e., throughfall, evapotranspiration, soil infiltration, vertical movement of water in the soil, and the accumulation of surface and subsurface runoff) [44]. Accordingly, Noah LSM provides these details as input for the routing modules. Ref. [45] describes a solution for spatially weighted sub-grid disaggregation–aggregation to coordinate the LSM with the terrain routing grids in WRF-Hydro [46].

2.3. Main Parameters Used in Hydrological Models

WRF-Hydro provides various hydrological parameters related to soil characteristics and watershed runoff, thereby controlling the total water volume and the temporal distribution of river discharge. Adjusting such parameters can yield potential improvements in the simulated data, depending on the study region. Some of the most important parameters in this study are denoted as infiltration factor (REFKDT), surface retention depth (RETDEPRT), surface roughness (OVROUGHT), and channel roughness (MannN) [47]. These four empirical parameters are dimensionless, often without an interpretable physical meaning, and are suggested to be adjusted through model calibration [48].

REFKDT is a parameter corresponding to saturated hydraulic conductivity for sandy clay soil and can significantly influence surface infiltration and the partitioning of total runoff into surface and subsurface flow. An increase in REFKDT leads to a reduction in surface runoff [48]. The equation involving REFKDT is defined by Equation (2):

$$K_{dtref}=K_{dt}\times \frac{K_{ref}}{K_{sat}}$$

(2)

where K_dtref is the parameter REFKDT, K_sat is the saturated hydraulic conductivity, K_ref is the saturated hydraulic conductivity for sandy clay soil, and e K_dt is a constant for calculating the maximum soil infiltration rate (I_max) as derived from Equation (3):

$$I_{max}=\frac{P\times \left\{\frac{d_{tot}\times \left[1-e^{\left(-K_{dt}\times ∆t\right)}\right]}{P+d_{tot}\times \left[1-e^{\left(-K_{dt}\times ∆t\right)}\right]}\right\}}{∆t}$$

(3)

where P is the effective precipitation intensity, d_tot is the total soil water depth (m), and t is the time interval.

The surface retention parameter (RETDEPRT), controlled by a surface retention depth scale parameter (RETDEPRTFAC), serves a function similar to that of REFKDT, indicating the amount of water the soil can store. The elevation level on each grid segment is an amalgamation of native infiltration surplus, the quantity of water streaming onto the grid segment from superficial flow, and the discharge from subterranean water flow. The quantity of surface height above the retention depth is accumulated as input to the stream channel and is effectively “discharged” to the channel routing routine. Enhancements in the RETDEPRT scaling coefficient in the channel cells can promote more vigorous local infiltration adjacent to the river channel, leading to more saturated soils that more accurately mimic riparian conditions.

The parameter for overland flow roughness (OVROUGHRT), controlled by the scaling parameter OVROUGHRTFAC, is determined by land use type and affects the downstream velocity of overland flows. It is a parameter used to represent how rough or smooth the land surface is, which affects the speed and efficiency at which water flows over it in response to precipitation. A higher value of OVROUGHRT makes the surface rougher, causing water to flow more slowly over it. Since hydrological models are often calibrated for a specific region or watershed rather than a specific vegetation type, it is recommended to modify the OVROUGHRTFAC scaling factor [13].

The MannN parameter reflects the influence of channel roughness on discharge, defined as Manning’s roughness coefficient. It is based on Manning’s equation, which is one of the most commonly used for calculating flow over free surfaces, as is the case in this WRF-Hydro study conducted in watershed channels. This parameter is used as the resistance formulation and requires the specification of the terrestrial flow roughness parameter [13,17].

In this study, a sensitivity test is conducted with the infiltration factor (REFKDT) while keeping the default values for the other parameters. After identifying the optimal value for the infiltration factor, tests for the surface retention depth scaling parameter (RETDEPRTFAC) begin, using the best REFKDT value obtained in the previous test and default values for the other two parameters. Similarly, tests are performed for the scaling parameter (OVROUGHRTFAC) and then for the MannN parameter [19,20,44,47,49].

2.4. Statistical Evaluation Metrics

The model’s performance in simulating discharge is evaluated based on statistical measures commonly used in other studies with WRF-Hydro. The indices employed include the Nash-Sutcliffe efficiency (NSE), root mean square error observations ratio (RSR), where RMSE is the root mean square error, percent bias (PBIAS), and the correlation coefficient (R) [15,20,50,51,52]. These statistical indices are defined by Equations (4), (5), (6), and (7), respectively, given by:

$$NSE=1-\left[\frac{{\sum }_{i=1}^N{\left(O_i-M_i\right)}^2}{{\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2}\right]$$

(4)

$$RSR=\frac{RMSE}{{STDEV}_{obs}}=\frac{\sqrt{{\sum }_{i=1}^N{\left(O_i-M_i\right)}^2}}{\sqrt{{\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2}}$$

(5)

$$PBIAS=\frac{{\sum }_{i=1}^N\left(M_i-O_i\right)}{{\sum }_{i=1}^N{O}_i}\times 100$$

(6)

$$R=\frac{\left\{{\sum }_{i=1}^N\left(M_i-\overline{M}\right)\left(O_i-\overline{O}\right)\right\}}{{\left\{{\sum }_{i=1}^N{\left(M_i-\overline{M}\right)}^2{\left(O_i-\overline{O}\right)}^2\right\}}^{\frac{1}{2}}}$$

(7)

where the indices $O_i$ e $M_i$ refer to the observed and model-simulated quantities for each day, respectively; $\overline{O}$ indicates the mean daily observed value over the entire period; ${STDEV}_{obs}$ is the standard deviation of the observed data; and N is the number of days in the simulation.

The NSE is a widely used performance index for assessing the quality of hydrological models. This parameter indicates how well the plot of observed data versus simulated data aligns. It ranges from −∞ to 1. The closer to 1, the better the model’s performance; values between 0 and 1 are generally considered acceptable levels of performance, and if the value is less than 0, the performance is deemed unsatisfactory [52].

The RSR is calculated as the ratio between one of the most commonly used error index statistics, o RMSE, and the standard deviation of the observed data (${STDEV}_{obs}$). The RMSE measures the differences between predicted and actual values, while the standard deviation measures the variability of observations around their mean. Thus, RSR provides a standardized measure of the model’s performance independent of the predicted variable’s scale. RSR ranges from zero, considered the perfect value, to a large positive value. The lower the RSR, the smaller the RMSE, and the better the model’s simulation performance [52].

The PBIAS index indicates the average tendency of the simulated data to be larger or smaller than their observed counterparts. The ideal value for PBIAS is 0; positive values indicate an average overestimation by the model, while negative values indicate underestimation [15].

The R value assesses the linear relationship between simulations and observations, ranging from −1 to 1. A value of −1 indicates a perfect negative correlation (when one variable increases, the other decreases in the same proportion), 1 indicates a perfect positive correlation (when one variable increases, the other also increases in the same proportion), and 0 shows no correlation [53].

3. Results and Discussions

3.1. Sensitivity Tests of Hydrological Parameters in Itajacá

The calibration procedure for hydrological parameters was based on several studies in the literature [16,17,19,20,49]. Sensitivity tests were initially conducted for the Manuel Alves Pequeno River using discharge data from the Itacajá station, focusing first on the REFKDT parameter that controls the hydrograph volume. Its default value is 3, with a feasible range from 0.1 to 10 [20]. This parameter can significantly affect the simulated discharge values, depending on climatic conditions and runoff generation mechanisms in different watersheds. Even in a physics-based model like WRF-Hydro, adjusting REFKDT and the other three empirical parameters is inevitable, a result of our inability to model all processes [48]. Based on the reviewed literature in this study, the values chosen for testing were 0.1, 0.3, 0.4, 0.5, 0.6, 1, 2, 3, and 4. Figure 4 illustrates the observed and WRF-Hydro-simulated discharges at the Itacajá station location for the specified REFKDT values from 13 March to 1 June 2018.

A significant impact of the infiltration factor (REFKDT) on model discharge results is observed, highlighting the sensitivity of this parameter, as demonstrated in [19] for the case in the Ouémé River, in [20] for the Upper Tana River Basin, both in West Africa, and by [47] on the Korean Peninsula. [49] conducted a calibration process similar to that of WRF-Hydro simulations. They used various values for REFKDT, REFDK, and OVROUGHRTFAC to identify the most effective parameter values for their basin, located in western Turkey. Their findings highlighted that REFKDT is the parameter with the highest sensitivity when it comes to regulating runoff reactions during rainfall events. They also observed that higher REFKDT values lead to a reduction in discharge. In the WRF-Hydro model, the water entering the channel, affecting the discharge, is primarily derived from the amount of precipitation exceeding the infiltration capacity at each grid cell [17]. As illustrated in studies in the Mediterranean region [16] and in the Daqinghe Basin in northern China [17], it is generally observed that an increase in REFKDT results in a decrease in surface runoff and, consequently, a reduction in simulated discharge. This trend is visible in the graph for the period from 1 April to 15 April 2018, where REFKDT = 0.1 shows a peak discharge close to 1400 m³ s⁻¹, REFKDT = 0.5, around 600 m³ s⁻¹, and for REFKDT values of 2, 3, and 4, the peak is around 200 m³ s⁻¹.

Table 2. Statistical indices of the comparison between observed discharge at the Itacajá station and the simulated discharge with REFKDT sensitivity tests.

REFKDT	0.1	0.3	0.4	0.5	0.6	1	2	3	4
NSE	−2.30	0.44	0.03	0.45	0.47	0.46	0.47	0.41	0.36
RSR	1.82	0.75	0.99	0.74	0.73	0.70	0.73	0.77	0.80
R	0.67	0.74	0.74	0.74	0.78	0.73	0.69	0.65	0.61
PBIAS	50.43	6.94	6.94	10.03	17.82	−7.31	−6.68	−11.58	−12.56

presents the statistical coefficients NSE, RSR, R, and PBIAS when comparing simulated and observed discharge for each REFKDT value.

This finding aligns with [54], who emphasizes the WRF-Hydro model’s capability to adeptly manage and simulate conditions of peak flows. Incorporating these advanced features of the WRF-Hydro model, our work strategically portrays a detailed and enriched understanding of flood peaks and their dynamic behaviors. Such precise modeling, supported by robust empirical evidence, amplifies the reliability of our hydrological forecasts and analyses, making them especially pertinent in navigating and managing scenarios characterized by extreme hydrological events.

The PBIAS index generally indicates a decrease in average discharge with the increase in REFKDT, representing maximum overestimation (PBIAS = 50.43) for REFKDT = 0.1 and the highest magnitude of underestimation (PBIAS = −12.56) when REFKDT = 4. The REFKDT values of 0.5, 0.6, 1, and 2 yielded similar statistical indices. REFKDT = 0.6 was selected as it exhibited the best correlation (R = 0.78), one of the highest NSE values, and the lowest error (RSR = 0.73). The results of the statistical metrics are also interesting for REFKDT = 0.3, as it showed a lower PBIAS and reasonable values for other metrics. However, we did not use this value because, during the second peak flow, around 22 March 2018, the simulated discharge showed a significant underestimation.

With the selection of REFKDT = 0.6, values of 0.1, 0.5, 1, and 3 were analyzed for the parameter RETDEPRTFAC. Its default value is 1, and it has a function similar to REFKDT, potentially influencing the calculation of surface runoff. This is because RETDEPRTFAC serves as a factor of maximum retention depth; before that, the flow is routed as overland flow [13]. Figure 5 presents the observed discharge values at the hydrological station and the simulated values for the mentioned RETDEPRTFAC parameters.

The values of RETDEPRTFAC 0.1, 0.5, and 1 evaluated in [17] showed virtually no influence on their discharge results. In this study, there is a slight variation in the simulated discharge at these values, specifically at peak points around 10 April 2018, suggesting that the influence of the RETDEPRTFAC parameter on channel discharge is not evident.

Table 3. Statistical indices comparing observed discharges at Itacajá station with simulated discharges from sensitivity tests of RETDEPRTFAC.

RETDEPRTFAC	0.1	0.5	1	3
NSE	0.45	0.38	0.47	0.46
RSR	0.67	0.78	0.73	0.62
R	0.75	0.69	0.78	0.79
PBIAS	13.64	12.79	15.82	−2.60

presents the statistical coefficients derived from the data represented in the plot.

In the WRF-Hydro modeling system, it is generally expected that surface runoff, and consequently, discharge should decrease with an increase in RETDEPRTFAC. This trend was not consistently observed except for the highest value, RETDEPRTFAC = 3, which displayed the lowest PBIAS index (−2.60), indicating the only underestimation of the simulated data relative to the observed data.

As observed in [17], or in cases in the Kaidu River Basin in China [55], in the Eastern Black Sea and Mediterranean regions [18], where discharge results showed nearly no response to RETDEPRTFAC, the statistical indices demonstrate lower variability compared to the other parameters, suggesting that RETDEPRTFAC has weaker sensitivity in simulating discharge among the parameters tested in this study. This may be due to significant variations in altitude and terrain in the study area, as suggested in [17], where little accumulation occurs on steep surfaces. Despite its limited influence, the default value, RETDEPRTFAC = 1, yielded the best correlation (R) and an NSE closest to 1. Figure 6 presents the discharge rates observed at the Itacajá station and those simulated by WRF-Hydro with sensitivity tests for OVROUGHRTFAC.

Just as in the study conducted in the Upper Trinity River Basin, located in North-Central Texas, USA [48], infiltration processes were more sensitive to REFKDT and OVROUGHRTFAC. Consequently, significantly altering the discharge data at certain OVROUGHRTFAC values changes the shape of the hydrographs in these cases, showing peak discharges at separate times and values. At the maximum point of the observed simulations, for example, around 10 April 2018, the discharge is greater when OVROUGHRTFAC is 0.1 than when this parameter is 0.5, which in turn is greater when OVROUGHRTFAC = 0.8.

Table 4. Statistical indices comparing observed discharges at the Itacajá station with simulated discharges from sensitivity tests of OVROUGHRTFAC.

OVROUGHRTFAC	0.1	0.2	0.4	0.5	0.6	0.8	1
NSE	0.44	0.62	0.57	0.63	0.69	0.61	0.47
RSR	0.75	0.61	0.66	0.61	0.56	0.62	0.73
R	0.77	0.79	0.80	0.80	0.83	0.79	0.78
PBIAS	13.26	10.35	10.06	10.38	4.99	6.00	17.82

shows the statistical coefficients for all parameters.

It is observed, as in [17], that this parameter had a quite significant impact on improving the statistical indices when varying OVROUGHRTFAC relative to its default value, except when OVROUGHRTFAC = 0.1. In these cases of improvement, the NSE, which is the statistical index used in all WRF-Hydro flow studies, improved by a scale of 21% to 47%. The other statistical indices also improved considerably, making it, along with REFKDT, the parameter that most influenced the flow simulation in this study.

Lastly, after analyzing the transmission of surface runoff to the river network, the transport of water along the channels was assessed through the sensitivity test of the channel roughness parameter (MannN). This transport can also affect the shape of the hydrograph. The channel geometric properties were set to their default values because there are no available channel cross-section data for the region concerning stream orders. In WRF-Hydro, channel properties such as average channel base width (Bw), initial water depth (HLINK), channel slope (Ch SSlp), and Manning’s coefficient (MannN) are introduced based on each stream order. The default values of the channel parameters are provided in

Table 5. Default values of channel parameters.

Stream Order	Bw	HLINK	Ch SSlp	MannN
1	1.5	0.02	3.0	0.55
2	3.0	0.02	1.0	0.35
3	5.0	0.02	0.5	0.15
4	10	0.03	0.18	0.10
5	20	0.03	0.05	0.07
6	40	0.03	0.05	0.05
7	60	0.03	0.05	0.04
8	70	0.10	0.05	0.03
9	80	0.30	0.05	0.02
10	100	0.30	0.05	0.01

.

The adjustment of the Manning coefficient is carried out by multiplying each initial MannN value by a scaling factor, thereby preserving the spatial patterns of the parameters during calibration. This scaling factor becomes the calibration parameter, generating MannN values for all the river channels. In this study, the values were established to cover a range between 0.2 and 2, including increments of 0.2 [19]. Figure 7 displays the simulated flow rates with their respective MannN parameters and the flow rate measured at the hydrological station.

When the MannN coefficients ranged from 0.2 to 0.8, the simulated flow rate peaks were significantly higher than the values obtained when MannN equals or exceeds 1. Flow rates exceeding 400 m³ s⁻¹ were not observed within this range, indicating a stabilization of the flow rate peaks.

Table 6. Statistical indices comparing observed flow rates at the Itacajá station and those simulated with sensitivity tests for MannN.

MannN	0.2	0.4	0.6	0.8	1	1.2	1.5	1.8	2
NSE	−0.81	0.10	−0.92	0.32	0.69	0.61	0.61	0.58	0.57
RSR	1.35	0.95	1.39	0.82	0.56	0.62	0.62	0.65	0.66
R	0.65	0.81	0.69	0.77	0.83	0.79	0.79	0.78	0.77
PBIAS	26.32	25.09	26.27	15.08	4.99	8.36	10.79	11.29	10.28

displays the statistical indices for each observed MannN value.

The more pronounced overestimation displayed graphically when MannN is less than 1 is corroborated by the PBIAS index, which shows higher results in these cases, reaching over 5 times the smallest value obtained when MannN is 1 (PBIAS = 4.99). All other statistical indices were better at the default value of the channel roughness parameter when MannN = 1.

Therefore, the best performance in flow rate simulation with WRF-Hydro in this study was achieved when REFKDT = 0.6, RETDEPRTFAC = 1.0, OVROUGHRTFAC = 0.6, and MannN = 1.0. The statistical indices NSE = 0.69, RSR = 0.56, R = 0.83, and PBIAS = 4.99 indicate superior results. According to [51], the range of NSE > 0.65 and RSR ≤ 0.60 is considered good for the evaluation of hydrological models, and the PBIAS index ≤ ±10% is considered exceptionally good. In their best flow rate simulation with WRF-Hydro in coupled mode, [19] obtained a correlation index R = 0.84, [44] achieved R = 0.83 and NSE = 0.7, and [20] reported NSE = 0.71 and RSR = 0.51.

3.2. Results at the Other Two Hydrological Stations

In the Jacaré and Goiatins stations, which measure the flow of the Vermelho and Manuel Alves Grande rivers, the same rain event from the Itacajá station was analyzed, covering the period from 13 March to 1 June 2018. Since it is the same simulation, the standard choice of hydrological model parameter values was maintained.

Table 7. Statistical indices comparing observed flows at the Jacaré and Goiatins stations with simulations using sensitivity tests for parameters REFKDT, RETDEPRTFAC, OVROUGHRTFAC, and MannN.

Jacaré Station					Goiatins Station
REFKDT	NSE	RSR	R	PBIAS	REFKDT	NSE	RSR	R	PBIAS
0.1	−0.61	1.27	0.64	24.54	0.1	0.06	0.97	0.75	8.87
0.3	0.25	0.87	0.61	−3.80	0.3	0.51	0.70	0.74	−12.08
0.4	0.17	0.91	0.61	−3.80	0.4	0.50	0.70	0.74	−12.13
0.5	0.33	0.82	0.59	−13.79	0.5	0.46	0.73	0.70	−12.67
0.6	0.37	0.79	0.63	−9.43	0.6	0.50	0.71	0.75	−13.96
1.0	0.29	0.84	0.57	−16.58	1.0	0.51	0.70	0.75	−16.08
2.0	0.17	0.91	0.48	−22.55	2.0	0.38	0.79	0.68	−21.29
3.0	0.17	0.91	0.48	−23.31	3.0	0.33	0.82	0.65	−22.39
4.0	0.18	0.91	0.51	−24.80	4.0	0.36	0.80	0.66	−20.67
RETDEPRTFAC	NSE	RSR	R	PBIAS	RETDEPRTFAC	NSE	RSR	R	PBIAS
0.1	0.30	0.83	0.60	−5.78	0.1	0.50	0.71	0.75	−10.18
0.5	0.19	0.90	0.50	−14.48	0.5	0.43	0.75	0.69	−17.17
1.0	0.37	0.79	0.63	−9.43	1.0	0.50	0.71	0.75	−13.96
3.0	0.38	0.78	0.64	−14.86	3.0	0.50	0.71	0.79	−17.44
OVROUGHRTFAC	NSE	RSR	R	PBIAS	OVROUGHRTFAC	NSE	RSR	R	PBIAS
0.1	0.16	0.91	0.59	1.93	0.1	0.50	0.71	0.73	−8.46
0.2	0.27	0.85	0.56	−16.95	0.2	0.50	0.70	0.73	−10.75
0.4	0.30	0.84	0.58	−13.79	0.4	0.49	0.71	0.73	−14.34
0.5	0.38	0.79	0.63	−9.85	0.5	0.56	0.66	0.77	−12.79
0.6	0.39	0.79	0.64	−9.60	0.6	0.61	0.63	0.80	−13.78
0.8	0.33	0.82	0.60	−14.75	0.8	0.49	0.71	0.71	−10.16
1.0	0.37	0.79	0.63	−9.43	1.0	0.50	0.71	0.75	−13.96
MannN	NSE	RSR	R	PBIAS	MannN	NSE	RSR	R	PBIAS
0.2	−0.91	1.38	0.42	7.10	0.2	−0.06	1.03	0.61	−4.65
0.4	−0.12	1.06	0.54	−4.86	0.4	0.05	0.98	0.56	−13.87
0.6	−0.33	1.15	0.38	−8.13	0.6	0.54	0.68	0.78	−18.37
0.8	0.36	0.80	0.62	−11.92	0.8	0.61	0.62	0.82	−16.39
1.0	0.39	0.79	0.64	−9.60	1.0	0.61	0.63	0.80	−13.78
1.2	0.31	0.83	0.58	−14.14	1.2	0.45	0.74	0.70	−9.44
1.5	0.30	0.84	0.57	−12.96	1.5	0.43	0.75	0.68	−7.90
1.8	0.32	0.82	0.59	−12.63	1.8	0.44	0.75	0.68	−7.66
2.0	0.26	0.86	0.54	−11.49	2.0	0.43	0.75	0.67	−7.39

shows the statistical parameters for the results of these simulations.

Similarly, in the analysis conducted at the Itacajá station, the same values assigned to REFKDT, RETDEPRTFAC, OVROUGHRTFAC, and MannN converged to the best statistical indices in most cases of the simulations that were analyzed with the Jacaré and Goiatins stations. Furthermore, the improvement in results was more pronounced when changing the parameters REFKDT and OVROUGHRTFAC.

The statistical coefficient PBIAS is often lower as the hydrological model parameters are larger, especially for the REFKDT parameter, indicating a decrease in surface runoff and a consequent reduction in simulated flow. Figure 8 shows the simulations compared with the Jacaré and Goiatins stations.

Indeed, it is noticeable that lower REFKDT values exhibit higher peak points compared to higher values, which show flow results close to each other. It can be observed in the first two peaks of the observed flow. The simulations, in the majority of cases, have their values underestimated. However, at the last peak, the simulated flows are overestimated.

Despite the three stations belonging to the same watershed and the same simulation tests being carried out, it is evident that the statistical indices of the simulations compared to the Jacaré and Goiatins stations are lower than those of the simulation conducted with the Itacajá station. This could be related to an increased disparity in the distribution of rainfall in both spatial and temporal dimensions, which worsens the WRF-Hydro simulations. This could result from errors in the directional data, in this case, the GDAS-FNL data used in the WRF, which indirectly affects the WRF-Hydro’s performance in flow simulation [17].

3.3. Error in Simulated Precipitation

The values of simulated precipitation, when compared to the values of measured precipitation using the MERGE product in the drainage network region of the Rio Manuel Alves Pequeno (where the Itacajá station is located), exhibit a more uniform distribution of rainfall over time than in the drainage network regions of the Rio Manuel Alves Grande and Vermelho (where the Goiatins and Jacaré stations are located, respectively). This directly affects the difference between the peaks of simulated and observed flows and the timing at which they occur.

Figure 9 displays the daily averages of precipitation obtained from the MERGE product within the drainage network area of each river, as well as those simulated by the WRF within the same area. It also shows the observed flows at river gauge stations and simulated flows in the WRF-Hydro model under the best-performing conditions, characterized by REFKDT = 0.6, RETDEPRTFAC = 1.0, OVROUGHRTFAC = 0.6, and MannN = 1.0.

It can be observed that the simulated precipitation in the observed data presented in Figure 9b,c is underestimated, while Figure 9a exhibits a seemingly smaller overestimation. It is also noteworthy that the inherent uncertainty in climate predictions resulting from the WRF simulation, conducted with GDAS-FNL input data, has impacted the results of the WRF-Hydro simulation.

Table 8. The average values of accumulated simulated and observed precipitation across all grid points within the respective drainage network of each river, along with the relative errors.

River	Observed (mm)	Simulated (mm)	ER
Manuel Alves Pequeno	517.55	572.73	10.66
Vermelho	441.45	304.56	−31.01
Manuel Alves Grande	407.55	290.26	−28.78

displays the relative errors of accumulated simulated and observed rainfall events within the drainage networks of the three analyzed rivers.

Indeed, the drainage network area of the Rio Manuel Alves Pequeno demonstrates the best results in precipitation simulation among the other two areas, with the smallest absolute value of ER = 10.66%. Conversely, the WRF model exhibits the poorest performance in simulating rainfall accumulation for the Rio Vermelho, with the highest relative error. This discrepancy could potentially impact the results of flow simulation in WRF-Hydro. Therefore, it is essential to emphasize that the calibration of WRF-Hydro parameters in this study cannot eliminate the errors in rainfall simulation originating from input data, directly influencing flow results. However, even with errors stemming from precipitation and considering the possibility of improving data through further study of WRF-Hydro in Brazilian regions, the strong correlation observed between simulated and observed flows when modifying hydrological parameters can serve as a starting point for using this hydrometeorological model in future research related to water resources management in Brazil.

3.4. Implications for Water Resources Management

Our study contributes valuable insights into the hydrological dynamics of the Manuel Alves Pequeno, Vermelho, and Manuel Alves Grande rivers, providing a foundation for improved water resources management practices in the MATOPIBA region. The accurate representation of hydrological parameters, such as soil infiltration, surface retention depth, and land surface roughness, is critical for making informed decisions regarding water allocation, especially in areas experiencing climate change and rapid urban and agricultural expansions.

Water resource management is a crucial aspect of ensuring sustainable water allocation, addressing the challenges posed by climate change, and expanding urban and agricultural areas. Accurate representation of hydrological parameters by numerical models such as WRF-Hydro is essential for making informed decisions regarding water resources [13]. By accurately simulating streamflow, the research serves as a benchmark for future hydrological modeling efforts in the MATOPIBA region, guiding water resource management strategies to meet the demands of a changing climate and evolving land use patterns. The efficient simulation of streamflow is a prerequisite for any effective water management plan. For example, [56] assessed a water security threat index for the western region of Bahia, an area within the MATOPIBA region and one of Brazil’s primary areas for agricultural exploitation. The authors estimated that 66% of agricultural enterprises effectively manage the water resources in the region. However, their results suggest that threats to water availability in the region will significantly increase by 2040.

Environmental exploitation, stemming from Brazil’s international significance as a food exporter, climate change, land use/land cover (LULC) changes, and desertification processes are key factors exacerbating this issue. LULC and climate change can significantly affect surface runoff, which is a key component of water resource management. Models that investigate the changes in hydrological pathways caused by climate and human activities can provide valuable frameworks for understanding and managing surface runoff [57,58]. Analyzing the impact of climate change on net irrigation water demand can help assess the future water resources in areas equipped for irrigation [59]. Therefore, considering the use of modern hydrometeorological models like WRF-Hydro can assist decision-makers, whether in the public or private sectors, in better managing the available water resources.

Climate change significantly impacts water resources and agriculture, highlighting the crucial need to understand how climate shifts interact with human activities in assessing water availability and food production [60]. In the bigger picture, global climate change poses threats to farming by messing with natural resources, especially rainfall patterns. Extreme downpours linked to climate change boost rainfall erosivity, causing more soil erosion. On the other hand, less rain creates water shortages for farming. Realizing how vital it is to model climate conditions under global change, particularly where rainfall data are on the rise, a study focusing on the MATOPIBA region becomes a big deal.

Ref. [61] projected rainfall erosivity in the Tocantins-Araguaia River basin for the 21st century under IPCC AR5 scenarios (RCP4.5 and RCP8.5). Using downscaling from four CMIP5 global climate models via the Eta regional climate model, the study showed a drop in erosivity for both scenarios due to less rainfall. But, even with this drop, the erosivity values are still significant, shouting out the need for robust soil conservation practices. Plus, the lower rainfall hints at less water for long-cycle crops and more uneven and less intense rainfall. In this complicated mix of climate challenges, our results with WRF-Hydro rainfall and streamflow simulation in the MATOPIBA region stand out, giving vital insights for managing water resources and planning agriculture to mitigate future scarcity.

Coupled hydrologic models, such as WRF-Hydro, are valuable tools for estimating water balance, analyzing groundwater levels and surface water discharges, and extrapolating to different scenarios and large basins [62,63]. In hydrological studies, the availability of data is vital. Data on rainfall, runoff, and other hydrological parameters are essential for accurate modeling and understanding of the dynamics of river basins [64]. These models can provide detailed information for water resource management and help analyze the impacts of meteorological and LULC changes.

Additionally, our research serves as a benchmark for future hydrological modeling efforts in the MATOPIBA region, guiding water resources management strategies to meet the demands of a changing climate and evolving land use patterns. While the manuscript may not explicitly delve into the policy and management aspects, the accurate simulation of streamflow is a prerequisite for any effective water management plan.

4. Summary and Conclusions

This study investigated the simulation of streamflow from 13 March to 1 June 2018, using the WRF-Hydro model in three rivers within the Tocantins/Araguaia sub-basin, between the Sono and Araguaia rivers in the MATOPIBA region. The analysis focused on assessing the influence of hydrological parameters, including the infiltration factor (REFKDT), surface retention depth (RETDEPRT), surface roughness (OVROUGHRT), and channel roughness (MannN). Sensitivity tests of these parameters were conducted concerning discharge data from the Itacajá station, located on the Manuel Alves Pequeno River. The parameter values selected for the Jacaré and Goiatins stations, which measure the flow of the Vermelho and Manuel Alves Grande rivers, respectively, followed the same pattern.

The study found that REFKDT and OVROUGHRTFAC were the parameters that had the most significant impact on the statistical indices obtained from the comparison between observed and simulated discharge. In many cases, the statistical coefficient PBIAS was lower when the model hydrological parameters were larger, particularly in the case of the REFKDT parameter. This suggests a decrease in surface runoff and a consequent reduction in simulated discharge.

The best discharge simulation performance was achieved when REFKDT = 0.6, RETDEPRTFAC = 1.0, OVROUGHRTFAC = 0.6, and MannN = 1.0. The highest statistical indices obtained were NSE = 0.69, RSR = 0.56, R = 0.83, and PBIAS = 4.99 in the simulation compared to observations at the Itacajá station. At the Jacaré station location, the statistical indices for discharge simulation were lower due to the less uniform temporal distribution of rainfall compared to the regions within the drainage network of the other two hydrometric stations.

These results indicate that the fully coupled WRF-Hydro modeling system demonstrates a good tendency in river discharge simulation, highlighting the considerable influence of model hydrological parameters and temporal rainfall distribution at each observation location. Future works involve improving results by analyzing additional model variables and incorporating other resources, such as artificial intelligence. Similar parameterizations will also be applied to different periods of intense rainfall in this region, aiming to effectively assess water resources in MATOPIBA basins and other Brazilian regions.