Abstract
Quantitative precipitation forecasts (QPFs) were obtained from ensembles of the weather and research forecasting (WRF) model for the Iguaçu river watershed (IRW) in southern Brazil. Thirty-two rainfall events between 2005 and 2010 were simulated with ten configurations of WRF. These rainfall events range from local to synoptic scale convection and caused a significant increase in the level of the Iguaçu river. In the average, the ensembles yielded up to 20% better skill than single WRF forecasts for the events analyzed. WRF ensembles also allow estimating the predictability through the dispersion of the forecasts providing relevant information for decision-making. Phase errors of ensemble forecasts are larger than amplitude errors. More complex microphysics parameterizations yielded better QPFs with smaller phase errors. QPFs were fed to IRW hydrological model with similar phase and amplitude errors. It is suggested that lagged QPFs might reduce phase errors.
1. Introduction
Hourly QPF is still a challenge in atmospheric sciences and weather forecasting. Uncertainties in numerical weather prediction (NWP) are caused by incomplete and inaccurate observations when assimilated and amplified by model’s dynamical and physical approximations. They can be reduced by ensemble forecasts studied since the 60s with statistic of NWP outputs [1]. Different initial conditions and physics are used at the same starting time to produce distinct solutions of the true future state of the atmosphere. The goal is to improve forecasting and provide uncertainties of each forecast [2].
NWP ensembles use several member alternatives, mass and energy perturbation, lagged analysis, data assimilation, and physical parameterizations. ECMWF operational ensemble system [3] uses a combination of forecast members with physical and dynamical perturbations to spread the range of solutions and to minimize ensemble mean errors. Physical solutions are provided by different microphysical and cumulus parameterization while dynamical solutions are carried out wind and temperature perturbation fields with the singular vector method. It estimates the growth of small perturbations over a period of time [4]. NCEP short-range ensemble forecast [5] was implemented operationally with a 10-member ETA in 2001. It was updated in 2009 [6] with 3 ETA members with BMJ cumulus scheme, 3 with Kain-Fritsch cumulus scheme, 3 with RSM/SAS and Ferrier microphysics, 2 with RSM/RAS and Zhao microphysics, 5 with WRF-NMM, and 5 with WRF-ARW. Eight of the WRF members use initial condition (IC) from GEFS (global ensemble forecast system) to increase the spread growth rate [7]. Environment Canada uses 20 members for both global and regional ensembles with Kalman’s filtering for dynamical perturbations and a set of physical parameterizations. In this research, similar RMSE performance was obtained with improved reliability in the regional one by more spreading in forecasts [8]. Four operational global numerical weather prediction models, ECMWF, JMA, NCEP, and GFS and UKMO, are used simultaneously with appropriate weights to obtain the multimodel ensemble (MME) technique on real time at India Meteorological Department, New Delhi [9]. For the Indian domain studies, giving weights to different models improves the QPF up to five days and bias-correct technique can produce better results [10].
At operational centers of meteorology, usually the forecasters use two or more models. They often analyze the model skill at the end of the forecast time lapse to verify whether or not phenomena such as cyclones, cold fronts, or convective systems are reproduced. This method can be used to select important meteorological features leading to severe weather. They might be useful to produce perturbations with better severe weather forecast skills than static perturbations [11].
QPFs depend upon appropriate microphysical and cumulus parameterizations and other larger scale phenomena. Models should be able to specify the water content by convection of individual cells with intense updrafts and downdrafts and stratiform regions composed of older cells with vertical velocities less than 1 m s−1. Furthermore, microphysical schemes account for a number of hydrometeor types resulting from condensation, evaporation, freezing, melting, deposition, and sublimation [12–14]. A myriad of physical parameterizations goes into forecast models to reduce errors and to improve QPF skills. WRF version 3 [13] has 14 microphysical schemes and 7 cumulus parameterizations. Multiphysic members can improve forecast skills of severe weather introducing [14].
Precipitation is the most important input variable to hydrological models. QPF for small and medium watersheds is still a major challenge, but recent studies indicate improvement on flood forecasting [15]. In general, the forecasts are more skillful for moderate precipitation amounts than those for either light or heavy precipitation [16]. Hydrological models have simpler physics than atmospheric ones and can run much faster. So, if you have computer capabilities for running an atmospheric model ensemble, it is most likely that hydrological multiple parallel runs using precipitation simulations are feasible [17]. Results of the rainfall-runoff model driven by the QPFs calibrated with methodology based on analogs methods revealed a beneficial impact on the reduction of missed events for the autumn and spring seasons of the years 2003–2008 [18]. Accordingly [19], the bias can be reduced by recalibrating the intensities of the ensemble mean using the probability density function (PDF) of rainfall values of each ensemble member. Even week-2 forecast, an improvement of about 18% on average is obtained for both streamflow and reservoir inflow forecasts using ensemble QPFs [20].
The Iguaçu river watershed (IRW) could be modelled with fuzzy, neural network, stochastic, and distributed models with compatible streamflow forecasts skills [21]. Usually, these IRW models run on measured precipitation and streamflow only. This paper shows results of the QPF ensembles to improve streamflow forecasting. One-hour resolution ensemble QPFs are used for this purpose for small and medium size watersheds.
2. Methods
2.1. Iguaçu River Watershed
The IRW is within Paraná state, Brazil (Figure 1(a)). This river stream is 1320 km long and flows from Curitiba city to Foz do Iguaçu. In its sloped channel 10 reservoirs were constructed for hydroelectric power generation. The IRW total area is 67,000 km2 and the area of the study is its first 24,000 km2, in which no reservoirs were constructed (Figure 1(b)) and wherefore the river has a natural flow. Altitudes vary from 943 m upstream, 765 m at the outlet of the study area, to 164 m downstream IRW where Itaipu hydroelectric power plant was constructed. It is the most important hydroelectric power plant in southern hemisphere, with Iguaçu falls downstream from it, one of the new seven wonders of nature.
(a)
(b)
QPF and streamflow simulations were carried out upstream from the IRW reservoirs. The watershed head is located at a large urban area of about 1,200 km2 and continues downstream on farmlands and forest. IRW terrain is quite irregular so the time of concentration of surface runoffs ranges from 6 hrs to 24 hrs. Several flood episodes occurred in the last thirty years related to ENSO (1983, 1992, 1997, 2005, 2007, 2010, and 2011) with major social and economic impacts. Each flood episode has specific consequences to urban and rural planning. Monitoring and forecasting of river stage remains an important activity for the protection of life and property. Noteworthy, the upstream IRW ends at União da Vitória city (UVC) where streamflow varies from 200 m3 s−1 to 400 m3 s−1 in winter and is above 1,000 m3 s−1 in summer. UVC is flooded with discharges greater than 1,300 m3 s−1 [22].
2.2. Weather Radar and Satellite QPE
Hourly rainfall measurements were obtained from a rain gauge network and estimates from an S-band weather radar within IRW (Figure 1(b)). The average distance between rain gauges in IRW is 27 km with higher density in metropolitan area (Figure 1(a)) of 10 km. The rain gauge network sensors are tipping buckets with resolution of 0.2 mm and inaccuracy of 0.5% for rainfall rates below 30 mm h−1 and 2% for above 120 mm h−1. The estimates of precipitation were obtained from radar using a traditional ZR relationship: where = effective radar reflectivity (mm6 m−3), = rainfall rate (mm hr−1), and , are adjustable parameters.
The classical Marshal-Palmer ZR relationship and [23] was used in this work. Both radar estimates and rain gauge-measured rainfall accumulation were integrated using the statistical objective analysis scheme (SOAS) developed by Pereira Filho et al. [24]. In this objective analysis scheme, the main goal is to minimize the analysis error by means of statistics obtained from rainfall measurements (gauges) and estimates (radar). Satellite rainfall estimates were obtained with IR channel data from GOES 13 and the CST technique [25]. The final rainfall field is a weighted average of gauge-measured, radar-estimated, and satellite-estimated rainfall accumulation [26]: where = integrated radar-gauge rainfall analysis at grid point ; = , the difference between radar-gauge and radar-only rainfall at grid point ; = integrated satellite-gauge rainfall analysis (mm) at grid point ; , the difference between satellite-gauge and satellite-only rainfall at grid point .
2.3. Ensemble Quantitative Precipitation Forecast (EQPF)
Hourly EQPF was obtained with WRF model, version ARW 3.1, applied in research and operation modes [27]. It integrates compressible and nonhydrostatic atmospheric momentum and cloud thermodynamics. Topography effects are prescribed by a terrain-following coordinate system.
A matrix of 10 WRF variants with a combination of physical schemes (Table 1) was run for 32 rainfall episodes between 2005 and 2010, representing a variety of mesoscale and synoptic conditions such as squall lines, cold fronts, and diurnal cycle convection over Iguaçu watershed. WRF was set to 42 vertical levels given by a tangent hyperbolic distribution with more levels close to the ground and tropopause. YSU planetary boundary layer, Dudhia, RRTM for short wave and long wave parameterization, Noah land surface model, MM5 similarity surface layer, and initial conditions specified by the 0.5-degree GFS global model were used in all WRF runs.
Recent results indicate that, with similar computer resources and model physics, the model resolution played a more important role than ensemble size when the forecast lead time was less than 4 days [28], even with coarse topography [29]. High resolution forecasts provide more detailed presentations of convective activity, but there appears to be little, if any, forecast skill on the scales where the added details emerge [30]. For this study the WRF was set to a 9 km horizontal grid resolution consistent with a good skill for mesoscale and operational computation resources.
WRF was integrated up to 48 h with a time step of 60 s and outputs every hour. The hourly average precipitation was input to the TopModel on a regular Cartesian grid model for each subwatershed.
2.4. TopModel
Hydrological simulations were carried out with the TopModel, a topography-based hydrological model [31]. It was selected based on its better treatment of watersheds with rapidly varying topography. The TopModel predicts the spatial distribution of water content as a function of the spatial variability estimated with an index of hydrological similarity and the mean overall water storage based on water balance at each time step [32].
TopModel assumes that the dynamic of the water cycle table is specified by a uniform subsurface runoff production per unit area over an area draining through a point and the hydraulic gradient of the saturated zone is given by the local surface topographic slope [31, 32]. The Iguaçu watershed was divided into nine subwatersheds (Figure 1) with outlets at river sections with streamflow measurements. Uniform soil and vegetation within each subwatershed are assumed. The topographic index () is given by where = cell of the subwatershed; = drainage area per unit length; = transmissivity of soil at the saturated zone; = slope angle of the surface.
The distribution of transmissivity with depth is an exponential function of storage. Deficit flow path is parallel to the landscape surface and the hydraulic gradient of saturated zone can be approximated by surface slope. The dynamics of the saturated zone can be approximated by successive steady state representations. According to Darcy’s law, the lateral flow at unit can be expressed as where = transmissivity of soil at the saturated zone, = slope angle of the surface, = storage deficit, = parameter indicating the decline of the transmissivity with soil, and = spatially uniform recharge rate to the water table.
The recharge from unsaturated zone to saturated zone for the whole catchment Qv is suggested as [31] where = unsaturated zone storage at unit ; total unit number for whole catchment; drainage area per unit length.
According to the storage changes in the unsaturated zone, the runoff owing to the saturation excess can be calculated. The total runoff is sum of runoff owing to the saturation excess and subsurface flow. The potential evapotranspiration was calculated by Penman method.
The slope angle was obtained from USGS GTOPO30 elevation data [33] at 1 km resolution.
The calibration of TopModel (Figure 2) was obtained with an independent data set between March 2003 and January 2004 when Iguaçu river stage level varied between 2000 m3 s−1 (May 2003) and 300 m3 s−1 (June 2003) [34].
2.5. Verification
The Brier skill score (BSS) was used to analyze probabilities of forecasts computing the average squared deviation between forecast probabilities and observational outcomes. BSS [35] can be divided into three terms: reliability, resolution, and uncertainty: where = number of forecasts in the category ; = probability of the class; = climatological base rate for the event occurrence; = observed frequency of the class; total of the forecasts issued.
The reliability term indicates how close the forecast probabilities are to the true probabilities. It is zero for a perfect forecast. The resolution term indicates how much probabilities differ from the climatology. It is no zero if forecasts differ from the climatologic average, so higher values are better. The uncertainty term depends only on the variability of the observations. It is one when the event occurs 50% of the time and is zero if the event always occurs [36]. For rare and frequent events, the uncertainty term is small. But, if the probabilities are at 50% of the time, the uncertainty is large.
QPF errors were estimated with the mean square error (MSE). It can be divided up into amplitude and phase errors or dissipation and dispersion errors, respectively [37]: where = observed precipitation observed; = forecasted precipitation; = standard deviation operator; = correlation coefficient between observed and forecasted precipitation.
Errors were estimated at one-hour time interval for all rainfall episodes.
3. Results
The EQPF WRF linked with TopModel was used to simulate 32 rainfall episodes over Iguaçu watershed in southern Brazil between 2005 and 2010. Most rainfall episodes started as deep convection followed by stratiform rainfall associated with cold fronts and squall lines. Cold fronts are more frequent in winter season with more even rainfall spatial distribution while squall lines and ordinary convection are in spring and summer seasons.
Mesoscale convective systems (MCS) are not common in the Iguaçu watershed region but can yield more than 100 mm day−1. Floods are in general caused by two or more consecutive precipitation events, prefrontal squall lines followed by cold fronts. The second rainfall shafts of the cold fronts yield uniform rainfall and tend to saturated soils and surface runoff quickly increases. Monthly precipitation average varies between 110 mm (winter) and 240 mm (spring) with an absolute maximum of 460 mm (October) and minimum of 3 mm (winter). Noteworthy is the range of monthly precipitation over Iguaçu watershed from 3.6 mm (August 2012) to 414 mm (August 2011).
Current hourly deterministic QPF has poor skill resulting in inaccurate precipitation forecasts input to hydrologic models. A byproduct of the ensemble can increase the skill of the forecasts and the set of members of WRF simulations allows uncertainty estimation, very useful operationally.
3.1. Ensemble Prediction System Diagrams (Epsgrams)
Ensemble prediction system diagrams (Epsgrams) are useful tools to estimate uncertainty of the precipitation forecasts. Epsgrams can easily extract minimum and maximum values, quartile tendencies, and convergence from the ensemble data.
The results indicated that the stratiform rainfall episodes tend to have smaller spread than convective ones though phase and amplitude errors persist. For instance, the ensemble forecast of a stratiform episode of September 11, 2005, the spread at = 20 h, was two times the observed precipitation average over the watershed (Figure 3(a)). In this case, the EQPF WRF can get successfully the phase and the amplitude of the average rainfall observed in the watershed.
(a)
(b)
(c)
(d)
Large forecast spread might be useful to detect the beginning of the rainfall event. The hourly precipitation forecasts (Figure 4(a)) for the August 01, 2009, prefrontal squall line (Figure 4(b)) ranged from 0 to 5 mm average over the watershed. But it was possible to detect the inset of precipitation. Some members predict zero mm and others do not detect it. The detection of this initial rainfall was important because next day after the squall line (Figure 4(b)) passed through the watershed and increased the river level, a cold front advanced over the watershed and produced stratiform rainfall (Figure 4(c)) over all Iguaçu watershed producing an intense flooding. EQPFs anticipate the rainfall stating time in a few hours (Figure 4(a)) but without delays afterwards which yield good hydrological performance for all time intervals.
(a)
(b)
(c)
Dispersion errors can be difficult to guidance. For instance, the April 25/26, 2010, rainfall episode had large uncertainty. An eastward bound squall line passed over the south of the watershed in the first day (Figures 5(b)–5(d)). The members predicted 0 mm hr−1 to 10 mm hr−1 with large scattering probably due to the position of the precipitation system. In the second day, a line of convection with trailing stratiform rainfall along the cold front occurred with maximum precipitation average over the watershed in good agreement with observations (Figures 5(e)–5(j)). But some forecast members yielded zero precipitation over the watershed. The members only agree at the end of the event.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.2. Brier Skill Score
BSS of hourly QPFs over the Iguaçu watershed of all rainfall episodes was below 0.25 with peaks between 0.3 and 0.45 (Figure 6), indicating that the forecasts were useful for all lead times. The resolution term did not capture the variety of events, indicating poor distinction compared to the climatology, probably because it used hourly comparisons. BSS maxima indicate forecast errors associated with phase errors given the hourly analysis is performed. The uncertainty is prevailing for all forecast times but in general with good scores for the first and the second day of forecast. They are also similar to others in the literature [38–40].
3.3. Deterministic and Probabilistic QPFs
Probabilistic forecasting is better in general than deterministic one in an operational use because it includes uncertainty, member’s convergence/divergence, ensemble means, and dispersion. A comparison between deterministic and ensemble mean simulations for Iguaçu watershed is shown in Figure 7. Probabilistic forecast yielded 10% to 20% better results. Ensemble means had a flat 0.8 hit score while deterministic ones varied from 0.5 to 0.8. At some times, the hit score of deterministic forecasts is equal to ensemble, never the opposite. Similar result was found for the southeast of the US using multimodel ensemble [41]. We assign this result to the better forecast of the phase of precipitation events (see Section 3.4).
3.4. Phase and Amplitude Errors
MSE, amplitude , and phase errors for deterministic and ensemble forecasts of all 32 rainfall episodes can be seen in Figure 8. The amplitude error of ensemble forecast was greater than deterministic in three periods of forecast 7–11 hr, 16–19 hr, and 28–32 hr (Figure 8). At the end of the first day and in the afternoon of the second day, the amplitude error of the deterministic forecast is much greater than ensemble mean one, probably due to the errors in simulating convection. Phase errors appear consistently greater than amplitude errors for both ensemble and deterministic forecasts. The use of ensemble mean could decrease the phase error for all times of forecast up to 48 hr. Phase errors cannot be reduced, but dispersion can indicate the inset of convection. It is useful information for hydrological simulation, especially for events of great impact [42]. Lagged runs reduced phase errors in some events, but it is not a rule. Furthermore, other members such as those used WSM6 and Morrison microphysics parameterizations also reduced phased errors in the same way as lagged members.
The WSM6 includes the prognostic of water substance variables: water vapor (), cloud water (), cloud ice (), snow (), rain (), and graupel (). Morrison is a two-moment scheme that predicts the same as WSM6 plus the number concentration (, ) and mixing ratios (, ) of cloud droplets and cloud ice . The ensemble members using WSM6 and Morrison produced deep convection faster than others providing more realistic hydrograph.
3.5. Precipitation Accumulation Statistics
The 24 hr precipitation accumulation bias of 32 events shows a large dispersion between the members and the QPE (Figure 9). These experiments reveal that the members tend to get the same bias signal, but, in some events, the range of dispersions is large both positive and negative. For the second day, the dispersion is greater than for the first day though the ensemble mean is improved. In many cases, the maximum and the minimum are similar in module. So, this dispersion is an evidence of differences among microphysical schemes and the importance of ensemble forecasts rendering uncertainty for numerical predictions.
(a)
(b)
The root-mean-square error (RMSE) and correlation (CORR) for 24 hr and 48 hr precipitation on Iguaçu watershed (Table 2) show important differences between members. Morrison, Lin, and WSM6 were the best microphysical schemes while Kain-Fritsch (KF) and Grell 3D cumulus schemes reached better skills than Betts-Miller-Janjic and Grell-Devenyi. The skill of the ensemble mean has been performed up to individual members carrying out accurate forecasts according to hourly simulations.
3.6. Hydrological Simulations
Hourly QPFs were input to TopModel to forecast streamflow at Iguaçu river for all 32 rainfall episodes. Precipitation scenarios were important to anticipate the increase in river stage. Iguaçu watershed has a time of concentration between 6 hrs and 24 hrs so the precipitation that falls on the basin could lead to a good prediction of the river’s level with a well-calibrated hydrological model.
The TopModel runs outputs using ensemble QPFs’ scenarios for January to May 2007 are shown in Figure 10. At the beginning of increase streamflow, the WRF member with WSM6/KF physics yielded the best phase in stage level increase. Both the ensemble mean and the deterministic forecast delay it. This result is an important guidance to deepen the analysis of the impact of the forecast and can serve as warning to go to further simulations. Comparing the calibration (Figure 2) with the example on Figure 10 we may conclude that the errors in the simulation of hydrological model come, largely, from the weather forecast and then improvement of the QPF skill has a direct impact on the streamflow prediction. The reducing bias 10–20% founded in the ensemble mean experiments will contribute to increase the predictability of the river flow.
Amplitude and phase errors of streamflow are shown in Figure 11 for both deterministic and ensemble precipitation inputs. Phase errors are greater than amplitude errors. This indicates that TopModel simulations are affected by both precipitation and streamflow measurements (memory), but the catchments can respond quicker than precipitation inputs. This difficulty is amplified in downstream points during extreme events, like hurricane [43].
4. Conclusions
An ensemble WRF with different physics and lagged analysis was used to QPF for a middle size watershed in southern Brazil. The ensemble forecast was compared to deterministic ones yielding better results for all 32 rainfall episodes. Ensemble means are better than deterministic forecasts with an overall skill 10% to 20% higher.
It also allowed estimating predictability, very useful to operational applications. Depending on the dispersion of forecasts, one can further simulate excluding members with poor skill or feed in more recent observations and analysis.
The results indicate that deterministic forecasts can be replaced by ensemble ones with several advantages to forecasters. Results with complex microphysics as WSM6 with graupel prognostic and Morrison scheme that include number concentration and mixing ratio of cloud droplets and ice prognostics improved the skill of the ensemble mean probably by the better solution of convection.
This work indicates that forecasts are dominated by phase errors for all rainfall episodes that also affect hydrological simulations with TopModel. Lagged runs might reduce phase error, but it is not a rule. Members with complete microphysics such as WSM6 and Morrison obtained better skill in precipitation forecasts, including reducing phase errors.
The TopModel simulations were affected by both precipitation and streamflow measurements (memory), but the catchments can respond quicker to the precipitation inputs. The EQPF members provided a large dispersion in the peak streamflow simulations and the extensive experiments show that the delay in the streamflow forecasts can only be reduced using more accurate QPF scenarios. So, the phase error in streamflow forecasts is a QPF dependent variable. The skill of the EQPF decreases drastically after the peak flow when a single member can solve the basic flow of the river.
As a suggestion to future work, members with data assimilation of physics characteristics of the atmosphere (radar data, e.g., [44]) and high-resolution cloud-scale simulations could be used to improve the skill of the ensemble.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to SIMEPAR and CNPq for the financial support Grants nos. 201341/2009-3 and 304033/2011-1 and to SIMEPAR for providing data sets.