TOPMODEL Hydrometeorological Modeling with Rain Gauge Data Integrated by High-Resolution Satellite Estimates. A Case Study in Muriaé River Basin, Brazil ()
1. Introduction
Floods are a natural phenomenon, but when occurring in inhabited areas they become natural disasters. Brazil is among the most affected countries, with 10,444 events recorded between 1991 and 2010, which resulted in more than 38 million people affected [1]. Urban expansion associated with socioeconomic problems acts in the increase and intensification of natural disasters related to floods. This is due to changes in the original hydrological cycle, such as increased soil waterproofing, deforestation, erosion, and structural interventions in rivers [2] [3]. The population (mainly low income) increase in hazardous areas enhances the magnitude of natural disaster’s impact [4]. The effects on the affected population can be diverse and quite significant: direct deaths, destruction of buildings, evictions, disruption of dikes and dams, and obstruction of roads [5]. The economic damage caused by natural disasters is also significant, as it is estimated that in 2008 alone it was approximately US$ 1 billion [6].
In Brazil, the South and Southeast regions, the latter where the Muriaé river basin (hereafter, MRB) is located, is the most affected by natural disasters caused by hydrological phenomena [6]. Severe floods were recorded in the MRB in 1997, 2008, and 2012 [7]. Most recently in 2020, intense precipitations resulted in severe floods and landslides, which 84 cities located at or near the MRB the state of emergency was decreed with 7 deaths and 40,000 unsheltered people [8].
One alternative to better understand hydrological phenomena is through hydrological models. They make it possible to study the impact of land-use changes in a watershed and the prediction of flood events [9]. In a study in the Paraíbado Sul basin (in which the MRB is located) Caluan & Cardoso (2020) [10] used a hydrological model to conclude that an increase in forest cover would decrease the frequency of very high flow events in the rainy season and also decrease the frequency of very low flow events in the dry months. Hydrological models can be defined as the representation of processes that occur in a watershed which allows the prevention of the consequences of the different occurrences concerning the observed values [2]. Beven (2001) [11] defines hydrological models as a means to estimate hydrological variables in space and time, as a support for decision-making related to water resources such as flood forecasting.
An example of a hydrological model is the Topography-based Hydrological model (TOPMODEL) [12], a model developed for small basins in humid regions of the United Kingdom whose main parameter is the topography of the basin. TOPMODEL is a model classified as semi-distributed and conceptual that aggregates the advantages of the complexity of a distributed model with the simplicity of concentrated and unique parameters for the entire basin [13]. It was successfully used in several hydrologic studies in small basins (i.e. <500 km2) [13] - [18] and medium and large size (e.g. 800 and 27,000 km2) [19] [20].
TOPMODEL’s primary input variable is precipitation. The calibration efficiency and the result of any hydrological model are dependent on the accuracy of the measurement of this variable and, therefore, should be estimated and measured with the slightest possible error [21]. However, obtaining reliable and representative data of precipitation measures can be an overly complex task due to spatial heterogeneity and the different temporal scales of rain events [22].
Rain gauges installed on the surface directly measure precipitation. However, they have limited spatial variability representation in regions with a low equipment density [23]. Besides, tipping bucket pluviometers may underestimate the measurement in situations of very intense precipitation [24] and strong winds [25].
Indirect measurements via remote sensing such as satellite estimates are an alternative. These methods allow for greater space-temporal resolution, but as an indirect measurement of precipitation, they contain several sources of errors and uncertainties [26]. Satellite quantitative precipitations estimates (SQPE) are not effective at estimating stratiform rainfall and fail to detect convective systems due to data sampling [23]. The Climate Prediction Center Morphing Technique (CMORPH) [27] is an example of an SQPE product that uses data from passive microwave (PMW) and infrared (IR) sensors onboard satellites to produce precipitation data with high spatial-temporal resolution.
One solution to reduce precipitation errors is to use objective analysis to integrate rain gauge data with remote and spatially distributed estimates of weather radars or satellite data [21]. The Statistical Objective Analysis Scheme (SOAS) [28] was developed to combine the advantages of each measurement system to produce a precipitation analysis field with minor errors and, therefore, improving the performance of hydrological simulations. Several studies have already used SOAS to integrate precipitation data and apply them to hydrological models with satisfactory results [17] [29] [30] [31] [32]. In Pereira Filho et al. (2018) [33], a sixth-order polynomial equation of CMORPH 24-hour precipitation correlation between two pixels as a function of the distance between them was adjusted. The correlation function was adapted for the CMORPH precipitation time series from 2000 to 2015 over the entire Brazil territory. This equation made it possible to apply the SOAS method in this study.
Since TOPMODEL was developed for small-sized basins with a temperate climate, a minority of its applications have been to medium-sized basins with a tropical climate. When applying TOPMODEL to a hectare scale tropical climate basin in French Guiana, Moliĉová et al. (1997) [34] concluded that the results could be susceptible to the temporal and spatial resolution of the precipitation data. In a 73 km2 basin in Southeast Brazil, Rocha Filho (2010) [17] observed that only three parameters showed high sensitivity with the result (m, lnTe, and td) and that the low sensitivity of the other parameters can mean a lack of adherence between the observed and calculated data [35]. When applying the model in a humid climate, 380 km2 basin in Nigeria, Campling et al. (2002) [16] changed the dynamics of the topographic index so that the water table was not parallel to the terrain. The calibration results were satisfactory and the subsurface flow was the most crucial process in generating the hydrograph. The authors also noted that intense convective precipitations were not properly represented by the network of rain gauges used, which corroborates the need to use precipitation measurement systems with high spatial resolution.
Because it is a simple model that each parameter (except topographic index) has a unique value for the whole basin (i.e., lumped model), TOPMODEL is a viable alternative in the operation of a flood warning system. Thus, this study aims to study the effectiveness of TOPMODEL for hydrological forecasting in the MRB, a medium-sized basin and tropical climate, a situation with a limited number of applications. For this purpose, SQPE from the CMORPH product and analyzed precipitation that integrated CMORPH data with rain gauge measurements with the SOAS method was used
2. Material and Methods
2.1. Study Area
The MRB is located in Southeast Brazil within the states of Minas Gerais (45% of the area) and Rio de Janeiro States (55%) with a drainage area of 8126 km2 (Figure 1). The Muriaé river is the last main tributary of the Paraíbado Sul river flowing to the Atlantic Ocean.
In this study, hydrological simulations were performed for Carangola and Patrocíniodo Muriaé subbasins (hereafter referred to as CLA and PMU respectively). Carangola stage level gage is at the main tributary of the eponymous river with a drainage area of 742 km2. The Patrocíniodo Muriaé stage level gage is
Figure 1. (A) Map of Brazil divided by federative units with an emphasis on the southeast region. (B) Southeast region of Brazil with an emphasis on the MRB location. (C) Map of the MRB (thick black line) with emphasis on CLA (light grey) and PMU (dark grey). Main rivers are represented by thin black lines. Stage gages are represented by a black triangle. Longitudes are represented on the abscissa axis and latitudes on the ordinate axis.
at the Muriaé river and has a catchment area of 2990 km2.
The Digital Elevation Model (DEM) used was the product of the Shuttle Radar Topography Mission (SRTM) operated by the National Aeronautics and Space Administration (NASA). The images were available at 30-meter horizontal spatial resolution and 1-meter vertical resolution. In Table 1 the land use and soil type of CLA and PMU are presented. The MRB terrain (Figure 2) is characterized by rugged terrains with topographic features of 1500 m altitude in the northwest quadrant and flat plains in the southeast quadrant.
The region climate is classified as hot and humid with a dry winter, Köppen class Aw [37]. Temperatures are high throughout the year, with a rainy summer (December to March) and a dry winter (June to August) when monthly rainfall averages are less than 50 mm (Figure 3). The average annual rainfalls are 1400 mm (PMU) and 1200 mm (CLA), while the average flow values are 48 m3∙s−1 (PMU) and 13 m3∙s−1 (CLA).
2.2. CMORPH
The Climate Prediction Center Morphing Technique – CMORPH [27] is an SQPE with a high spatial-temporal resolution (~8 km and 30 min). The method is based on precipitation rates estimated from PMW sensors. The PMW images are obtained from three sensors, located at polar orbit satellites which yield images for a specific location with a 3-hour temporal resolution (Table 2).
The IR images on board of geostationary satellites (METEOSAT and GOES) are used to increase the spatial-temporal resolution of the final product through wind vectors from the Cloud System Advection Vector (CSAV) method [41]. The wind vectors advect the PMW rainfall rate fields in space and time. Consecutive PMW images are advected forward and backward in time. In morphing, two images from the same time step are composed to generate the final product.
Table 1. Land use and soil type for CLA (second column) and PMU (third column). Sources: http://sigaceivap.org.br/siga-ceivap/map and https://www.ibge.gov.br.
*Classes defined by [36].
Figure 2. Digital Elevation Model (DEM) of the MRB obtained from the SRTM product. Carangola and Patrocínio do Muriaé stage gages are represented by black triangles. Altitudes are represented by the color chart (greyscale) in meters above sea level. Black lines represent the main rivers. Longitudes are represented on the abscissa axis and latitudes on the ordinate axis.
Figure 3. Climatology of Patrocínio do Muriaé (black) and Carangola (grey) hydrometeorological stations. Monthly mean precipitations (mm) are represented by vertical bars (downward on the left ordinate axis) while monthly mean discharges (m3 s−1) are represented by lines (upward on the right ordinate axis).
Table 2. Satellite sensors used for CMORPH. Sensors, satellites, microwave frequencies, and references on respective precipitation estimates are shown.
Overland, both temperature and emissivity are highly variable, which makes the microwave background signal unpredictable and makes it difficult to use images from emission-based frequencies for precipitation estimates [42]. Therefore, over continents, PMW-based precipitation estimates rely solely on the scattering of higher frequencies microwave from frozen hydrometeors.
2.3. SOAS
An alternative to minimize precipitation errors is through the SOAS method [21]. An efficient interpolation technique [43] was initially developed to integrate radar precipitation estimates and rain gauge measurements. Recently, the methodology was also used to integrate SQPE [33].
The analyzed precipitation at a given grid point is obtained by adding precipitation estimates by remote sensing (radar or satellite) at the point and the sum of the product of the weights of the differences between rain gauges measurements and those estimated by remote sensing [17]. In this case, the error variance will be the sum of the squares of the differences between the measured values and the actual values. Therefore, the weights of the analyzed values will be dependent on the error variance of the precipitation estimates. The analysis is performed so that the error variance is less than the error variance of the measurements. The SOAS equation [21] is:
(1)
where
is the analyzed precipitation at a point i of the grid [mm];
is the estimated precipitation by remote sensing at a point i of the grid [mm];
is the measured precipitation by the rain gauge at point k [mm];
is the estimated precipitation by remote sensing at point k [mm];
is the weight that will be set a posteriori; k is the number of rain gauges; (x, y) are the coordinates in UTM [km].
One premise of the method is that both rain gauge measurements and precipitation estimates are not correlated or biased. Thus, it is possible to obtain the weights through:
, (2)
where
is the correlation between the rain gauges k and i errors;
is the correlation between point i of the grid and the rain gauge k errors;
is the normalized observation error;
is the weight that will be set a posteriori.
The analysis error can be normalized by applying previously elaborated algorithms [44] in the error covariance matrices. The expected analysis error variance (NEXERVA) (
) can be expressed as:
(3)
NEXERVA is the ratio between the expected analyzed data error variance and the remote sensing estimates error variance. Thereby, it is possible to determine the spatial distribution of the errors in the analyzed field [21]. In general, NEXERVA is minimal at the rain gauge locations.
The matrix of the covariances of the estimated value errors used to normalize Equations (1) and (2) is the most crucial component of the SOAS [21]. This matrix has great importance in the accuracy of the analysis.
(4)
where
is the correlation between pixels k and l errors;
are the precipitation estimates by remote sensing in pixel k(l);
is the average long-term precipitation estimates by remote sensing in pixel k (l).
2.4. TOPMODEL
TOPMODEL is a rainfall-runoff model defined as deterministic, dynamic, physical-based parameters, and semi-distributed [12]. The model was developed for small humid basins of temperate climate [45]. The model is considered semi-distributed because only the parameter related to the topography, the topographic index, is spatially distributed. The topographic index (λ) can be expressed as λ = ln(α/tanβ), where α is the drainage area per unit of contour and tanβ is the terrain slope. High values of λ indicate points more likely to reach saturation and consequently more likely to generate runoff [45]. Points with equal weights of λ are treated as regions with similar hydrological responses.
Beven and Kirkby (1979) [12] and Beven et al. (1984) [45] described TOPMODEL as a physical-based conceptual model that represents surface and soil dynamics based on the relationship between discharge and storage established from flows in saturated zones due to slope. According to the authors, TOPMODEL combines the advantages of the simplicity of lumped models with the distributed effects of variable areas of contribution and flows through the drainage network. It is also possible to determine the parameter values from the knowledge of the physical characteristics of the hydrographic basin.
The transformation of precipitation into flow into a hydrological model can be divided into two stages: the soil water balance and the flow propagation through the basin [19]. In TOPMODEL, a series of three interconnected reservoirs represents the flow generation process and is the average response of the homogeneous capacity basin: the vegetation reservoir (S1), the soil unsaturated zone reservoir (S2), and the soil saturated zone reservoir (S3) [17]. For each time step, the reservoirs S1 and S2 are filled with the levels
and
respectively. The saturation deficit (S) of each time step is the difference between the maximum reservoir level S2 (
) and the instantaneous level (
).
This three reservoirs system generates four possible flows [19]: subsurface flow at the saturated zone (qs), vertical subsurface flow at the unsaturated zone (qv), surface runoff generated by infiltration excess (qh) (Horton mechanism), and surface runoff generated by saturation excess (qas) (Dunne mechanism) [46] [47].
In the next step, TOPMODEL propagates the surface runoff for the basin outlet hydrograph composition [48]. The main processes involved in the model are presented in Figure 4.
The model uses the premise that the simulation occurs after a dry period which implies the following initial conditions [19]:
· the unsaturated soil zone (S2) is completely dry (i.e.
);
· the discharge in the basin outlet is generated exclusively by the subsurface flow, represented by the parameter qs0.
The TOPMODEL’s parameters are presented in Table 3.
Figure 4. Schematics of TOPMODEL hydrological processes. S1: Vegetation reservoir;
: Maximum reservoir S1 level;
: instantaneous level of the reservoir S1; S2: unsaturated soil reservoir;
: The maximum level of the reservoir S2;
: instantaneous level of the reservoir S2; S: saturation deficit; S3: The reservoir of saturated soil zone; ETsup: Evapotranspiration that occurs on the surface; ETsoil: evapotranspiration that occurs in the soil; (1): precipitation; (2): effective precipitation; (3): saturated area runoff; (4): Horton mechanism runoff; (5): vertical subsurface flow; (6): surface runoff propagated with Clarke’s method; (7): subsurface flow.
Table 3. List of TOPMODEL’s parameters and their respective description.
2.5. Calibration with SCE-UA
For both subbasins, the same calibration methodology was used for all cases, the Shuffled Complex Evolution developed at the University of Arizona (SCE-UA) [49]. The methodology was developed to increase the probability that the best set of parameters found is a global minimum rather than a local minimum. The SCE-UA method was successfully used in the calibration of several hydrological models, such as NAM/MIKE 11 [50], Sacramento Soil Moisture Accounting (SAC-SMA) [51] [52], and TOPMODEL [53].
The probability of a global minimum was increased by the SCE-UA method. This method evolves multiple groups of parameter sets independently with the Competitive Complex Evolution (CCE) [49] [51]. After each round of the simulation, the parameter sets are rearranged to evolve with different partners. The simulation is repeated until the user-defined objective function of the best parameter set meets convergence criteria.
The SCE-UA calibration methodology was applied through the rtop library [54] of free software R (version 3.3.2). The objective function chosen was the Nash-Sutcliffe efficiency [55], which is given by:
, (5)
where NSE is the Nash-Sutcliffe efficiency;
is the simulated discharge at time step i;
is the observed discharge at time step i;
is the mean observed discharge of the whole period; n is the number of time steps.
For the simulations, a period of 4 years was selected between 2009 and 2013. To meet TOPMODEL’s premise that simulations must start in a dry period (i.e. absence of surface runoff), the month of September was defined as the beginning of the hydrological year. The hydrological years of 2009/2010 and 2010/2011 were chosen for calibration, while the hydrological years of 2011/2012 and 2012/2013 were selected for validation.
Daily rain gauge precipitation accumulation measurements between 1000 UTC (Coordinated Universal Time) of consecutive days and respective CMORPH precipitations estimates were used in the integration process with SOAS. Both rain gauges and discharge data were obtained from the Brazilian Hydrometeorological Network, which is operated by the Brazilian Water Agency (ANA). The SOAS method is described in detail in Pereira Filho et al. (2018) [33]. Evapotranspiration data were obtained from the closest weather station (Itaperuna, OMN #83695), operated by the Brazilian National Institute of Meteorology (INMET) in which evapotranspiration is measured with a lysimeter.
3. Results and Discussion
3.1. SOAS Applied to the MRB
Fifteen (15) rain gauges were selected (Figure 5(A)) for the SOAS method at the MRB. In the expected analysis error variance (NEXERVA) results (Figure 5(B)), it can be noted the influence of the rain gauges of the pluviometers, in which the greater the proximity to the observed data, the more significant the reduction of the error variance. The NEXERVA values for CLA (0.34) and PMU (0.38) subbasins show that the SOAS methodology reduced the CMORPH estimates uncertainties by approximately 65%.
3.2. Calibrated Parameters
The calibrated parameters of each simulation and the distribution of the topographic index of each subbasin are shown in Table 4 and Figure 6, respectively.
To comprehend the importance of each parameter, sensitivity tests were applied to all simulations (Figure 7 and Figure 8). The sensitivity analyses were
Figure 5. (A) Distribution of the fifteen (15) rain gauges used to apply the SOAS method for the MRB. Pluviometers are represented by crosses. Main rivers are represented by thin black lines. (B) NEXERVA values for the MRB. Longitudes are represented on the abscissa axis and latitudes on the ordinate axis.
Table 4. Calibrated TOPMODEL’s parameters for PMU and CLA with CMORPH and SOAS precipitation values.
Figure 6. Variation of the accumulated area (%) (y-axis) with the topographic index classes (x-axis) for PMU (black line) and CLA (grey line).
obtained through simulations in which only a given parameter has a variable value to determine its relevance in the calibration efficiency.
It can be noted that for all four simulations, three parameters presented the greatest sensitivity: qs0, m, and Srmax. The primary relevance of the parameter m coincides with previous studies [17] [45], which indicates the importance of the soil transmissivity variation with depth. The significance of the parameter qs0 suggests that the basin’s initial condition has great importance in the model’s performance. The high sensitivity of the parameter Srmax indicates the evapotranspiration importance. Higher values of Srmax for PMU are in agreement with the physical characteristics of its subbasin, which has greater forest coverage and less urban area (Table 1). The parameter vr (i.e. runoff speed) had higher values with SOAS, which indicates that in this system, precipitation tends
Figure 7. Sensitivity test of TOPMODEL parameters for PMU with CMORPH (blue) and SOAS (red). Parameter values are shown on the x-axis. Nash-Sutcliffe efficiencies are presented on the y-axis. Text: parameter – precipitation data system – subbasin.
Figure 8. Similar to Figure 7 except for CLA.
to be delayed with CMORPH, which makes greater speed necessary to compensate for the shorter time.
3.3. Simulations Performances
In PMU, simulations with CMORPH datasets (Figure 9) yielded similar performance for both calibration (NSE = 0.67) and validation (NSE = 0.66) (Table 5). Most peak flows were underestimated. Simulations with SOAS (Figure 10) presented higher performance for calibration (NSE = 0.60) in relation to validation (NSE = 0.39) (Table 5). Most peak flows were underestimated in the calibration and overestimated for validation for the hydrological year 2011/2012 and underestimated for 2012/2013.
In CLA, simulations with CMORPH (Figure 11) yielded a better result for the validation period (NSE = 0.61) over calibration (NSE = 0.54) (Table 5). In most cases, peak flows were underestimated. The fact that the validation efficiency was superior to that of calibration is an indication that the calibrated parameters satisfactorily represented the characteristics of the subbasin. In simulations with SOAS (Figure 12), the calibration period presented good efficiency (NSE = 0.85), with peak flows been well represented in the hydrological year 2010/2011 (Table 5). However, the simulations were not satisfactory for the validation period (NSE = −0.04). In the 2011/2012 hydrological year, the peak flows were overestimated approximately 100%, which explains the negative NSE value (i.e. the series mean would have been a better predictor).
For both subbasins, the highest NSE values were obtained with CMORPH data. However, in the analysis of high-level discharges (Table 6), simulations with SOAS resulted in a higher success rate (28% and 18% for PMU and CLA respectively) compared to CMORPH (18% and 13%).
Figure 9. Discharge (m3∙s−1) time series for PMU with CMORPH for observed (black line) and simulated (grey line) data. River stage level thresholds (attention, alert, and flood) are indicated by dashed lines. The date format is as follows: dd/mm/yyyy.
Figure 10. Similar to Figure 9 except for SOAS data.
Figure 11. Similar to Figure 9 except for CLA with CMORPH data.
Figure 12. Similar to Figure 9 except for CLA with SOAS data.
Table 5. NSE for TOPMODEL simulations of CLA and PMU with precipitation estimates from CMORPH and SOAS.
Table 6. Comparative analysis of TOPMODEL simulations for CLA and PMU with precipitation estimates from CMORPH and SOAS. First line: the total number of days in which the observed discharge was higher than the attention level discharge. Second line: number of days (percentage of days) in which the model underestimated the observed discharge. Third line: number of days (percentage of days) in which the model underestimated the observed had an accurate result. The fourth line: number of days (percentage of days) in which the model overestimated the observed discharge.
4. Conclusions
To better understand its application on medium-sized tropical climate basins, hydrologic simulations with TOPMODEL were performed on the Muriaé river basin with SQPE from CMORPH precipitation datasets and analyzed precipitation with SOAS that integrated rain gauges measurements and CMORPH estimates. Sensitivity analyses were performed to understand the relevance of each parameter to the results. For the overall period, simulated data with CMORPH had the best hydrologic simulation performance, while SOAS yielded the best simulation for high flow events (Table 6).
The results obtained in this study agree with other studies, in which only a small number of TOPMODEL’s parameters showed sensitivity. The parameter m relevance indicates the importance of understanding the soil transmissivity decay profile for tropical climate applications. Regarding evapotranspiration, the model showed physical significance. The more forested and least urbanized basin behaved so that the vegetation had a greater storage capacity (represented by the higher values of the Srmax parameter).
When compared to SOAS, simulated discharge values with CMORPH were significantly lower. Hypotheses for these findings are the inability of SQPE products to capture isolated convective systems due to sampling and the presence of precipitation produced by clouds without frozen hydrometeors, not detected by high-frequency PMW sensors (from 60 GHz to 300 GHz). For high stage level periods, the simulated flows with SOAS improved because of the recovery of the amplitudes of the precipitation field with the observed data applied to the high spatial resolution (8 km) precipitation estimates.
However, for the entire period, simulations with CMORPH resulted in better NSE values. That was an unexpected result since the SOAS method reduced the error variance of CMORPH estimates (Figure 5). Some factors that may have contributed to this result: short simulation period (four years); uncertainties in the rate curve (which are higher for extreme stage level values in which there are fewer measurements); possible lack of accuracy in some of the rain gauges. Since SOAS yielded notably worse results at the validation stage, it can indicate that the calibrated parameters were not representative of the subbasin and behaved as empirical values to better adjust to the input data.
Thus, TOPMODEL and these new high-resolution precipitation datasets can be used in MRB hydrometeorological forecasting, provided that the initialization of TOPMODEL starts in the dry season with high-resolution precipitation datasets. Therefore, this hydrometeorological forecast improves flood warning systems for civil protection.
Acknowledgements
The authors are grateful to CPRM for providing stage and rain gauge datasets. The authors would like to thank CPC/NCEP for providing high-resolution CMORPH datasets. They also would like to thank three anonymous reviewers who enriched this manuscript. The second author is supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under grant 302349/2017-6.