Coupled Hydrologic-Hydraulic Modeling of the Upper Paraguay River Basin

The conceptual model presented in this study has two main components: (1) simulation of the basin and part of the Paraguay River tributaries by means of a rainfall-runoff model with a simplified flow routing method; and (2) simulation of the main drainage network by a full 1D hydrodynamic model.

The large-scale, distributed hydrological model MGB-IPH was used as the first component of the conceptual model. It was developed for use in large South American basins with scarce data, and is fully described in Collischonn et al. (2007a). The drainage basin is represented by square-grid elements connected by channels, each one of these elements contains a limited number of distinct grouped response units (GRUs), i.e., areas with similar combination of soil and land cover (Kouwen et al. 1993), similarly to the hydrological response unit presented by Beven (2001). Each element has dimensions of $0.1\times 0.1°$ (approximately $10\times 10\mathrm{km}$) in general; however, it depends on the basin size and model spatial discretization adopted. The MGB-IPH model consists of modules for calculating soil water budget, evapotranspiration, flow propagation inside an element, and flow routing through the drainage network. A soil water budget is calculated for each GRU, and runoff generated from the different GRUs in the element are then summed and propagated to the stream network using three linear reservoirs: base flow, subsurface flow, and surface flow. Streamflow is routed through the river network using the Muskingum-Cunge method, considering different length and slope for each river reach representing the connection between two given elements, which are automatically derived from digital elevation model (DEM) processing (Paz and Collischonn 2007; Paz et al. 2006).

The MGB-IPH model has been typically applied with a daily time step as this is commonly the time discretization of available rainfall data, although smaller time steps could also be used. Precipitation data are interpolated to the center of each model grid element at each time step using the inverse-distance-squared method, whereas a nearest neighbor method is applied for other meteorological variables. The calibration of the MGB-IPH model is achieved by changing parameters values while maintaining relations between them and landuse. The multiobjective evolutionary algorithm MOCOM-UA (Yapo et al. 1998) is employed for MGB-IPH automatic calibration considering three objective functions: volume bias ($\Delta V$); Nash-Sutcliffe efficiency index for streamflow (NSS); and Nash-Sutcliffe efficiency index for the logarithms of streamflow (NSSLOG). Several applications of the MGB-IPH model for hydrological modeling and streamflow forecasting at large basins have been carried out with acceptable results (Allasia et al. 2006; Bravo et al. 2009; Collischonn et al. 2007b; Collischonn et al. 2005; Tucci et al. 2003).

The second component of the conceptual model is a 1D hydraulic-hydrodynamic model, which was used to perform river flow routing along the main drainage network of the Pantanal region. The well-known Hydrologic Engineering Center–River Analysis System (HEC-RAS) hydraulic model (USACE 2004) was used. It solves the full Saint Venant equations using an implicit Preismann four point scheme of finite differences (Cunge et al. 1980), and Manning roughness coefficients are used to represent the resistance to flow.

The MGB-IPH model was applied to the entire Upper Paraguay River Basin with a spatial discretization of square-grid elements of 0.1° resolution and a total of 5,195 grid elements, with surface area ranging from 114.5 to $120.2{\mathrm{km}}^2$. A daily time step was used as precipitation, streamflow, and meteorological data are available in this time interval.

The river drainage network represented in the hydrologic model was automatically derived from the SRTM-90m DEM (Fig. 3), including flow direction and accumulated drainage area for each grid element and length and slope of river reaches connecting model elements. However, even using a stream burning procedure for preprocessing the DEM (Saunders 1999; Turcotte et al. 2001), manual corrections were necessary in the drainage network owing to the extreme flatness of the Pantanal.

The scarce precipitation and meteorological data (Fig. 3) were interpolated to the center of each grid element of the hydrological model at each time step. Soils were characterized using data available from the RADAMBrasil survey (Ministério das Minas e Energia 1983), PCBAP project (MMA 1997), and soil map published by the Food and Agriculture Organization of the United Nations (FAO) (1974, 1988). The soil types of the UPRB were grouped into eight classes according to their occurrence area and water storage capacity. Land-use classification was obtained by analyzing several Landsat7 ETM+ images, aiming basically at identifying areas occupied by crops, pasture, or native vegetation and seasonally flooded areas. According to available soil and land cover maps, soil types and land-use classes were combined. This resulted in a number of combinations, which were then grouped in 10 GRUs; this restriction in the number of GRUs aims at reducing model parameters and also to avoid representation of irrelevant or nonrepresentative physical characteristics for the scale of this study.

Table 1 summarizes the inputs and outputs of the hydrologic model for its application to the UPRB.

The hydrodynamic model was applied for flow routing along the Paraguay River and its main tributaries, representing its vast shallow lakes and divergent and endorreic drainage networks through a combination of reaches and storage areas approach.

As described in Paz et al. (2010), input data for the hydraulic model was prepared in a consistent and georeferenced database through GIS-based automatic procedures developed for dealing with the large amount of data provided by several different sources, with distinct formats and horizontal and vertical datum. Detailed cross-section profiles of the main channel were combined with elevation values extracted from SRTM-90m for characterizing the floodplains, maintaining the link between hydraulic data and spatial location. Following these procedures, all the geometric data relative to the 24 river reaches, 12 junctions, 1,124 cross sections, and 11 storage areas were entered into the hydraulic model in a coherent and not time-expensive way.

Model calibration consisted in the setup of Manning’s coefficients ($n$). In this study, the coefficients calibrated by Paz et al. (2010) were used. These authors defined distinct values of $n$ for the main channel and floodplain for each river flow segment between two streamflow gauge stations by comparing observed and calculated hydrographs. This procedure was performed from upstream to downstream along the modeled drainage network manually varying $n$, but restricting them to the range reported in literature for these kinds of channels and floodplains. In total, there were 27 distinct river flow segments to have their corresponding $n$ coefficients calibrated. They run the HEC-RAS hydraulic model with a 12 hour-to-hour time step for the period of January 1, 1996 to December 31, 2000. Manning coefficients ranging from 0.02 to 0.035 were obtained for the main channel, and from 0.04 to 0.2 for the floodplain. Table 2 summarizes the inputs and outputs of the hydraulic model for its application to the UPRB.

A one-way and off-line coupling between the hydrologic and hydraulic models was adopted. The HEC-RAS hydraulic model fitted by Paz et al. (2010) was run with a 12-h time step for the same time period adopted by those authors (January 1, 1996 to December 31, 2000), but now using as input at the upstream boundary conditions the daily flows calculated by MGB-IPH model.

The entrance of each river into the Pantanal area was adopted as the upstream boundary condition in the hydraulic model (Fig. 4) as downstream from this point the flood wave is not correctly routed by the Muskingum-Cunge methodology used in the hydrologic model. Thus, the hydrographs calculated by the MGB-IPH model at these upstream boundary conditions were considered as the contribution of the Planalto region for the downstream river network modeled by the HEC-RAS model. In addition, streamflow generated in the Pantanal and Chaco areas through the rainfall-runoff transformation represented by the MGB-IPH model was considered as lateral inflows to the hydraulic model.

Model parameters were primarily set according with each GRU's characteristics, but some of them were adjusted by the automatic calibration process available for the MGB-IPH, restricting the parameter search between realistic physical limits. In spite of being a physical distributed model, as the number of distinct GRUs was set to a maximum of 10, because of computational limitations, some regions of the basin were not represented with enough details to correctly characterize the local hydrological processes. Consequently, slight differences among subbasins were allowed for varying the MGB-IPH model parameters related to GRU characteristics.

Different time periods were considered for calibrating each of the subbasins at the Planalto according to data availability. A total of 15 parameters of the MGB-IPH model was calibrated. The description of each parameter and its influence on distinct aspects of hydrological processes representation are described by Collischonn et al. (2007a). As a result of the calibration process, several Pareto optimal solutions were found for each gauging station, and a single solution was chosen among them, aiming at providing an acceptable tradeoff in fitting of the different parts of the hydrograph (Bastidas et al. 2002).

In general, the model fitted well as the NSS and NSSLOG coefficients were approximately 0.80 in most catchments in both calibration and validation periods (Table 3). The errors in volume between observed and calculated hydrographs were also acceptable, being close to 2–4% in most cases. Examples of the fit between calculated and observed hydrographs are presented in Fig. 5 for eight control points in the Planalto region, in which six of them represent the outlet sections of rivers flowing from Planalto to Pantanal.

The statistics presented previously should be understood in the context of data scarcity. The study used data from 86 streamflow and 92 rainfall gauges, which indicates a density of one discharge gauge every $2,953{\mathrm{km}}^2$ and one rainfall information every $2,760{\mathrm{km}}^2$ in the best and rare situation of having available data in a given day at all the stations simultaneously. These densities are too far from the World Meteorological Organization’s recommended ones for regions characterized by difficult data acquisition. Consequently, worst (best) results were achieved in the least (best) monitored basins. For example, it was not uncommon to have only one working rainfall station in the Aquidauana basin during most of the period, resulting in a density close to one rainfall information every $10,000{\mathrm{km}}^2$. In light of this, results obtained were considered acceptable.

The first interesting analysis of the results at Pantanal is the comparison to the results obtained in the former study of Paz et al. (2010), which used observed streamflow as an upstream boundary condition instead of coupling a distributed hydrological model. The statistics obtained after comparing observed and calculated hydrographs along the Paraguay River and tributaries were similar in both studies in terms of NSS and root-mean-square error (RMSE) (Fig. 6). Overall, the better the fit obtained by Paz et al. (2010), the smaller the difference between the statistics of those authors and this paper's. In other words, the difficulties in reproducing the observed flow regime at some streamflow gauge stations found by the discussed study, in which the Manning coefficients of the 1D hydrodynamic model were adjusted, were enlarged when the modeling of the rainfall-runoff processes over their contributing areas was included. This was expected because the uncertainties in precipitation estimates and in the rainfall-runoff modeling are now incorporated in the modeling. However, the similarity between the results achieved by this study relative to that obtained by the simplified approach of using observed discharges as upstream boundary conditions highlights the very reasonable results of the effort to model the entire UPRB.

The coupled model satisfactorily reproduced the observed flow regime at the tributaries, as illustrated by visual inspection of the hydrographs (Fig. 7) and according to the statistics obtained and showed in Table 4, in which MAE means mean absolute error and MRE means mean relative error. The NSS coefficients ranged between 0.7 and 0.9 for half of the available gauge stations and were less than 0.5 for only three of them (São João, São Jerônimo, and Porto Ciríaco). Along the tributaries, the MRE ranged between 12 and 22%, except at the Barão de Melgaço, São João, and Miranda stations, in which this measure was greater than 30%.

The major difficulty in reproducing flow regime along the Paraguay River tributaries is mostly because of uncertainties related to scarcity of available cross-section profiles. However, the proposed coupled hydrologic-hydraulic model was capable of reproducing the distinct observed flow regimes at the tributaries of the Paraguay River. At Barão de Melgaço, the Cuiabá River presents a marked seasonal flow regime, with the flood period occurring between October and May and peak flows reaching $1,600{\mathrm{m}}^3·\mathrm{s}$, whereas the recession flows are approximately $100{\mathrm{m}}^3·\mathrm{s}$. The São Lourenço River at the São José station also presents a marked seasonal flow regime, but this is characterized by smoother raising and falling limbs in comparison to the Cuiabá River at Barão de Melgaço, which presents more nervous oscillations in response to precipitation events. At São José, peak and recession flows are typically approximately 400 and $160{\mathrm{m}}^3·\mathrm{s}$, respectively, which mean flood peaks of just 1.5 times greater than recession flows. This ratio is 10 times greater for the Cuiabá River at Barão de Melgaço. However, downstream of this station, at Porto do Alegre, the flow regime of the Cuiabá River becomes very similar to that of the São Lourenço River at São José: smooth raising and falling limbs of flood hydrographs with relative small variation between flood peaks and recession flows. At São José Piquiri in the Piquiri River, a behavior slightly similar to this is observed but with more pronounced flood peaks. In turn, the flow regime at Miranda station in the Miranda River presents weak seasonality but with rapid and very pronounced response of flows to precipitation events. Lastly, at the Porto Ciríaco station in the Aquidauana River, the flow regime is very distinct to all those discussed previously, presenting a marked maximum value of $150{\mathrm{m}}^3·\mathrm{s}$ and slight seasonality. The observed hydrograph at this station presents only small peak flows, as during the major floods a huge volume of water spills over the main channel and inundates the floodplains, resulting in a slash of the raising limb at the value of $150{\mathrm{m}}^3·\mathrm{s}$.

All the complex flow regimes described were satisfactorily reproduced by the coupled hydrologic-hydraulic model, as indicated by visual comparison of observed and calculated hydrographs and by statistics shown in Table 4. Additionally, the remarkable changes in flow regime along each river flowing from Planalto to Pantanal were also reproduced by the proposed model. These changes are a consequence of the low conveyance capacity of river channels at Pantanal and the resultant floodplain inundation. For instance, compare the hydrographs between Rondonópolis (Fig. 5) and São José (Fig. 7); at the former, which represents the São Lourenço River outlet section from Planalto, seasonal peak flows range between 900 and $1,200{\mathrm{m}}^3·\mathrm{s}$, whereas at São José, the hydrograph is much smoother and peak flows are less than $450{\mathrm{m}}^3·\mathrm{s}$. Another notable example is along the 230-km reach of the Aquidauana River between the Aquidauana (Planalto) and Porto Ciríaco (Pantanal) stations, in which peak flows reduce from $700{\mathrm{m}}^3·\mathrm{s}$ to less than $150{\mathrm{m}}^3·\mathrm{s}$.

As obtained in the study of Paz et al. (2010), the worst results were found for the São João station, located in the Cuiabá River upstream the confluence of the São Lourenço and Piquiri Rivers. Along this river reach, the main channel capacity is reduced, and floods are diverted to the floodplain even at the smaller discharges; more than 50% of the annual volume is lost to the floodplain, and the peak discharge decreases from $2,000{\mathrm{m}}^3·{\mathrm{s}}^{-1}$ to less than $400{\mathrm{m}}^3·{\mathrm{s}}^{-1}$. Water in the floodplain may remain up to 3–5 months stored and evaporating before return in a diffused and not very well understood way through the complex system of secondary channels, shallow lakes, and ponds. However, the good model fit obtained at the Porto de Alegre station, for instance, $\mathrm{MRE}=18\%$ and $\mathrm{NSS}=0.73$, located 50 km downstream of São João suggests that the bias found at São João becomes less important after the inflow of the other tributaries of the Cuiabá River.

The results obtained at the Paraguay River show that the applied model was able to reproduce the marked seasonal flow regime and the typical relatively smoothed shape of hydrographs, as illustrated by Fig. 8. In Fig. 8, calculated and observed hydrographs at six stations along the Paraguay River are shown, depicting the changes in hydrographs according to flood routing from Cáceres up to Porto Murtinho, a reach of approximately 1,300 km. Along this flowpath, the Paraguay River receives significant contributions along its left margin, mostly from drainage basins of Cuiabá and Miranda Rivers. More importantly, these contributions occur by means of waters drained by both the main channel and floodplains, i.e.,  water flowing along floodplains of tributaries may contribute to channel flow along some reaches of the Paraguay River. On the contrary, along other reaches of this river, huge volumes of water spill over the main channel and inundate the floodplain. These multifaceted lateral water exchanges between channel and floodplain strongly make the hydrodynamic modeling of the Paraguay River difficult. Thus, the results obtained are very satisfactory as the applied model reproduced the variations on time and magnitude of flood peaks and recession flows along the Paraguay River, as indicated by statistics summarizing the comparison between observed and calculated hydrographs. At the Cáceres, Descalvados, Porto Conceição, Amolar, and Porto da Manga stations, NSS was greater than 0.74. The worst results were obtained at the São Francisco ($\mathrm{NSS}=0.50$) and Porto Murtinho ($\mathrm{NSS}=0.48$) stations. Along the Paraguay River, RMSE ranged between 52 and $510{\mathrm{m}}^3·\mathrm{s}$, and MAE ranged from 44 to $402{\mathrm{m}}^3·\mathrm{s}$. These values are apparently too large but correspond to MREs varying from 10 to 12%, except at São Francisco (21%) and Porto Murtinho (16%). The difficulty in reproducing observed flow regime three of them ( Francisco may be because of occurrence of contributions drained by Taquari River floodplains, which was not well-represented in the proposed model. At Porto Murtinho, the major reason for not achieving better results is the data scarcity for characterizing the contributions of the Bolivian part of the basin along the right margin of the Paraguay River. This was already suggested in the study of Paz et al. (2010), which indicated the secondary flood peaks in the observed hydrograph as being probably generated in the Bolivian part of Pantanal. On the contrary, an earlier study had concluded that the contribution of the Bolivian side of the basin was insignificant and could be disregarded (MMA 1997), considering the loss of water because of floodplain inundation and evaporation process and using a simplified hydrological balance procedure on the basis of available discharge data, i.e., data from streamflow gauging stations located in the Brazilian part of the basin. The study indicates, however, that for better fitting hydrograph volumes and timing along the Paraguay River downstream of Amolar, the contribution from the Bolivian part of the basin must be better quantified and represented in the model.